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PROGRESS REPORT NO. 8 

Proceedings of the Twenty-Fourth Seminar 

on 

Space Flight and Guidance Theory 


SUMMARY 

Progress reports of NASA sponsored studies in space 
flight and guidance theory are presented. The studies are 
made by several universities and industrial firms under 
contract to MSFC. This progress report reflects work done 
on the contracts during the period from April 1 , 1965 to 
December 31 * 1965 . The contracts are technically monitored 
by personnel of the Astrodynamics and Guidance Theory Division, 
Aero-Astrodynamics Laboratory, George C. Marshall Space Flight 
Center. 


INTRODUCTION 

This progress report contains eight papers whose subject 
matter lies within the areas of space flight and guidance 
theory. The papers have been written by investigators employed 
at universities and industrial firms under contract to MSFC. 

This report is the eighth of the "Progress Reports" and 
covers the period from April 1 , 1965 to December 31, 1965. 

The contributing agencies and their fields of major 
interest are : 


Stability of Dynamical Systems 

Trajectory Optimization 
Control Theory 


Brown University 
' General Precision Aerospace 
^ Grumman Aircraft 

f Southern Illinois University 
j The Boeing Company 
Auburn University 
r 

| Brown University 


The objective of this introduction is to review briefly 
the contributions of each agency. 

The first paper is concerned with finite time stability 
properties of periodic solutions of Hamiltonian systems. 

The contract for which it is a report of progress has for 
its principal objective the determination of that set of 
initial or injection conditions in which is included the 1 

initial conditions for a given (almost) periodic solution 
that will result in (almost) periodic solutions that lie 
within a prescribed "tube" of the given solution. A theoret- 
ical basis for this determination was given by Birkhoff. 

However, to obtain actual numerical estimates from the analyt- 
ical theory an enormous amount of algebraic manipulation is 
required even in the simplest problems. For this reason, a 
digital computer and an appropriate non-numeric computer 
language to perform the required manipulations were employed. 

As an initial application of the mechanization of the 
Birkhoff theory, the planar restricted three-body problem 
was chosen as a simplified dynamical model; the given 
solution in this model is a Lagrangian critical point. 

The second paper is concerned with discontinuous vector 
fields which are encountered in problems of feedback control. 

It begins with the observation that if X is a discontinuous 

vector field then the study of stability under perturbations 

e(t) is different if the perturbation enters the equation of 

motion as a summand in the argument of X, that is, as in the « 

equation 


x(t) = X(x(t ) + e (t) ) 

from what it would be if the perturbation were not a part 
of the argument, as in the equation 


x ( t ) = X(x ( t ) ) + e(t) . 


If X is continuous then this is not the case. Problems in 
feedback control lead to discontinuous vector fields in the 
form 


X(x) = F(x, m, ( x) ) 
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where p, is a control function. The author discusses the 
distinction between a classical solution of the equations 
of motion and a Filippov solution, a generalization of the 
definition of solution. He then shows as his main result 
that if a vector field X is stable with respect to measure- 
ment, then every classical solution is a Filippov solution. 

The third paper presents a survey of various approaches 
to the problem of estimating the domain of attraction of an 
equilibrium solution of a system of nonlinear autonomous 
differential equations. Based on observations resulting 
from this survey, the problem is reformulated as that of 
choosing optimally the Liapunov function from the space of 
positive definite quadratic forms. An estimate of the domain 
of attraction is then obtained as the solution of a minimi- 
zation problem. This approach to the problem has the advan- 
tages of being suitable for machine computation, of yielding 
estimates that are easily visualized and of being relatively 
insensitive to system dimension. Some preliminary numerical 
results are presented for the Duffing equation with damping. 

The fourth paper deals with the problem of obtaining a 
transformation technique which can be used to eliminate the 
control angles from the Euler-Lagrange equations to give a 
system of differential equations in the state variables and 
the Lagrange multipliers only. The problem arises in the 
study of trajectory optimization by classical calculus of 
variations techniques. In applying these techniques, certain 
Euler-Lagrange equations involving the control angles are 
encountered. In some cases these equations lead to a solution 
for the angles In terms of the Lagrange multipliers, and these 
solutions can be used to eliminate the control angles from the 
Euler-Lagrange equations resulting in a system of differential 
equations in generally desirable state variables and Lagrange 
multipliers only. The process, however, can be carried out 
more readily in some coordinate systems than in others. In 
this paper the technique for a general transformation of the 
state variables and their corresponding Lagrange multipliers 
from one coordinate system to another is discussed. The 
technique is then applied to a specific problem involving 
three-dimensional trajectory optimization. 

rp"h o -Pn -p-f-Tn o r>o *i <3 t*H +-Vi -Pn >0 ^ Cr’iCL 

the stability of the zero solution of the differential 
equation 

x^ + Pi(t)( n-1 ) + ... + P n-1 (t)x + p n (t)x - 0 
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which depend on the behavior of the real continuous functions 
p i (t), but not upon their derivatives. 

Recently, Ghizzetti obtained simple stability criteria 
for this problem. The particularly attractive aspect of 
these criteria is that they depend only on n constants which 
locate a family of hyperellipsoids in the n-dimensional space 
of the pj_(t). If the curve represented parametrically by the v 

Pi(t) is entirely contained within one of the hyperellipsoids, 
then the zero solution of the equation above is asymptotically 
stable. 

In this paper the author uses the second method of 
Liapunov to obtain stability criteria for the above equation 
which depend on only n parameters which determine a family 
of elliptic paraboloids in the n-dimensional space of the 
Pj_(t). It can be shown that these elliptic paraboloids 
completely contain the hyperellipsoids of Ghizzetti. A 
practical technique for the application of the stability 
criteria obtained is discussed and is applied to two examples. 

The objective of the sixth paper is to present a unified 
exposition of Liapunov's theory of stability that includes 
the classical Liapunov theorems on stability and instability 
as simple corollaries. The principal idea exploited in this 
paper was used by other investigators in the study of nonauton- 
omous functional differential equations. Of considerable 
importance is the possibility of extending these concepts to 
more general classes of dynamical systems, especially to some t 

types as defined by partial differential equations. 

A noteworthy contribution is Theorem 1 and its corollary. 

The theorem, which is concerned with the nonautonomous system 
x = f(t,x), explains precisely the nature of the information 
given by a Liapunov function; it shows that a Liapunov func- 
tion relative to a set G defines a set E which, under the 
conditions of the theorem, locates all positive limit sets 
of solutions x(t) of x = f(t,x) that for positive time remain 
in G. However, in order to use the theorem, there must be 
some means of determining which solutions remain in G. A 
corollary, a consequence of the theorem, gives one way of 
doing this and also provides, for nonautonomous systems, a 
method for estimating regions of attraction (domains of 
stability) . 

A limit set of Q is defined as the set approached by 
a solution x(t) of a system of differential equations as 
t-*oo. The points peQ are limit points. A limit set has an * 
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"invariance property" if all solutions x(t) which start at 
pen remain in fi as t-®°. It is pointed out that there are 
special classes of differential equations where the limit 
sets of solutions have, additionally, an invariance property 
and that this property permits a refinement and sharpening 
of Theorem 1, mentioned above, for these special classes. 

Because the paper is largely a survey of recent extensions 
of past investigations, formal proofs, except for corollary 6, 
are not given; but ample references and illustrative examples 
are provided for the reader who might wish to work out the 
proofs for himself. 

In the seventh paper an analytical solution of the 
Euler-Lagrange equations for the Lagrange multipliers for 
optimum coast trajectories is obtained. Similar solutions 
have been obtained by other investigators, but all of these 
solutions had singularities for orbits with zero eccentricity. 
The solution presented in this paper does not have such a 
singularity, but there is a numerical difficulty due to a 
removable singularity at unit eccentricity. An approximate 
solution, accurate near unit eccentricity, is given. This 
solution reduces to the exact parabolic solution for unit 
eccentricity. 
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Finite Time Stability of Periodic Solutions of 
Hamiltonian Systems 

by 

I. S. Bernstein 
Mathematics Department 


General Precision Aerospace 
Little Falls, New Jersey 


Abstract 


This study deals with finite time stability properties of periodic solutions of 
Hamiltonian systems. It attempts to answer questions such as what tube of 
initial conditions about the periodic solution will keep solutions in some pre- 
scribed tube about the periodic solution over some prescribed time interval. 

A theoretical basis for answering questions such as the one formulated above 
was given by Birkhoff [ 1 ] . However, to obtain actual numerical estimates 
from the analytical theory requires a great deal of algebraic manipulation 
even in the simplest of problems. For this reason, it was decided to employ 
a computer and an appropriate computer language to perform the required 
manipulations. 

As an initial application of the mechanization of the Birkhoff theory the 
system chosen was the planar restricted three-body problem and the solution 
chosen was a Lagrangian critical point. 

Basic Theory 

Let H(x,y) be the Hamiltonian of a dynamical system so that the equations 
of motion can be written as 

x = H (x,y) x = (x,,...,x ) 

v y ' 1,1 'In 

v = 1,...,n, (1) 

y v =- H x (*/y) y = (>V** ,/y n) 

V 

and let x^ = <P v (t), y^ = 0^(0 be a periodic solution of (1) of period 2ff. 

To study solutions in the neighborhood of this periodic motion we make the change 
of variables 

* = x - ip, y = y - 0. 
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Then (1) takes the form 


x v = H 7 (*' Y> 

7 V 

v-l,...,n, (2) 

7 V =-H~ (*/ 7/ 0 

v 


»w 

where H has period 2ff in t and 

H— <0,0,t) = hi_ (0, 0, t) — 0, v — l,...,n. 
'v v 


Thus the origin is a critical point solution of (2). The problem is therefore 
reduced to the study of solutions in the neighborhood of an equilibrium 
solution of a Hamiltonian system with an explicit periodic time dependence. 
A theorem of Birkhoff, [ 1 ] , is applicable to this problem. 


Theorem 1 Let the Hamiltonian H (x,y,t) of a dynamical system with 

an equilibrium point at the origin, be analytic in x and y, periodic in t 
of period 2ir, and thus representable in a convergent power series by 


H(x,y,t) = 



r 2'* 2n 


V 1 v 2 

, v „ (') x , x 2 ••• X - y 


v n V n+1 V n+2 


n *1 


^2 


2n 


V] + v 2 +...+ v 2n 2 


= H 2 (x,y,t) + Hg (x,y / t) +. . .+ H k (x,y,t) +. . . 

where H k (x,y,t) is a homogeneous polynomial in x,y of degree k with 
periodic coefficients of period 2 it. Let the 2n characteristic exponents, 

[2 ], associated with Hj be distinct and purely imaginary. As the system is 
Hamiltonian, they may be represented as, [3 ], 
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Xj,,,., \ ■ ^j/* •• ” x^. 

Furthermore let the exponents satisfy 

Assumption 1 m 1 X 1 +m 0 X.+...+ m X +m ,.^0 

II l 2 n n n + 1 ' 

for all integers m. such that 
n 

0< | m | = y l m 'l ^ N s 3 • 

— I 

1 = 1 

Then there exists a canonical change of variables 


x v f v U'^ 

v 1 , . . . ,n, 

y v = s v ( t' f ) 


( 3 ) 


where f^ and g^ are convergent power series without constant terms in 
the components of £ and tj with coefficients having period 2 V in t, 
such that the Hamiltonian in the new variables has the form 

h (£/ tj/ 1) = Hj (| 1 tj j ^ n n ) + h 2 U, n* 0/ 

where Hj is a polynomial with constant coefficients of degree N if N 

is even and degree N-l if N is odd in the variables z = £ r? , and 

v ^v 'v 

where Hj (£ , TJ / t) is a power series in £ y / T} y/ beginning with terms 
of degree N + 1 . With H in this form the Hamiltonian is said to be 
normalized up to order N. 
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Theorem 2 (Special Case) 


Let the Hamiltonian, Hk,y), of a dynamical system with an equilibrium 
point at the origin, be analytic in x and y and thus representable in a 
convergent power series by 


H(x,y) = 


«o 

5"= 

V r”- /V 2n 


V, V 0 V V , , V v 0 

12 n n+1 n+2 2n 


x j X2 • • • x 


n *1 ^2 “-yn 


y j+. . .+ Vjn 2 


H 2 (x,y) + Hg(x,y) +...+ H n (x,y) +..., 


where H n (x,y) is a homogeneous polynomial in x, y of degree n with 
constant coefficients. Let the 2n eigenvalues associated with H 2 be 
distinct, purely imaginary and represented by 


I n I n 


Furthermore let the eigenvalues satisfy 


Assumption 1 1 


m . X . + m-X 0 + ...,+mX *0 
II A i n n 


for all such integers m. such that 


0 < I m I 


= ^ | m. | ^ N s3. 


i=l 


Then the conclusion of Theorem 1 holds. Moreover there is no explicit time 
dependence in the change of variables (3) and thus also in the new Hamiltonian 
H (£, rj). 
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The usefulness of this theorem is the following. If we are studying solutions 
near the equilibrium point 4 = r) = 0^ then H- is of higher order than 

A. * 

and is discarded for the moment. The equations then have the form 


4 = (H.) (z) = f (H.) 

V S V 1 2 


t? v = -(H,) (z) = -t ?v (H 1 ) z 

9 V V 


v = l ;ilt/ n . 


If we multiply the first equation by rj^the second by 4 y j Q nd add, it 
follows that 


dt U v \) = 0/ v = l,...,n 


Thus, l y t? y = c^ (constant) so that (4) becomes integrable yielding 


£ = 4 (0) e' (H l*z t 

V V 


TJ = n (O)e’ 1 (, Vz (c)t 

V *V V 


v = 1 , . . . ,n . 


If we restrict ourselves to a large finite time interval and a suitable region in 
phase space it can be shown that the higher order terms previously truncated can 
be made small so that (5) is a close representation to the actual solution in this 
region. By use of (3) approximate solutions to the original problem may be 
obtained. For precise statements along these lines see [ 1] and [3], 

Rather than prove this general theorem we now illustrate how to carry out the norm- 
alization procedure for the particular Hamiltonian describing the planar restricted 
three-body problem and take for the periodic solution a Lagrange critical point. 
This problem falls under Theorem 2. 
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Application 


4 


The equations of motion for the planar restricted three-body problem in the 
rotating coordinate system are 


d^x 0 dy 2 
— _ - 2^ = a) x 

dt^ dt 


m 1 (x - x^ 

g 

r 1 


m 2 ( x - x 2 ) 

g / 

r 2 



+ 



2 

CO y 


mjy m 2 y 



where the gravitational constant has been set equal to one. 
If we set 

dx . , 
u = it ' y ' 


v = it + Wx ' 


then the Hamiltonian takes the form 


H -S 

(u 2 + V 2 ) + 

^ 1 "» 2 

w (uy - vx) . 

r 1 r 2 

(6) 

If we set |x^ 

-Xjl - d 

and 


d 

fn^ - «r»2 

. _/3d 


a 

m 1 + m 2 ' 

b — j" / 
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then the point 


x = a, y - b, u - - cob, v - t * J a / 

is an equilibrium point solution for the system. We introduce dimension- 
less variables T, q ]7 q 2 , Pj/ P 2 / in the neighborhood of this equilibrium 
point by 

T = wt, X = a + q^, y = b+q 2 d, u = - wb + p^d, v = u>a + p 2 <^d. 

The Hamiltonian (3) is now defined in a neighborhood of the equilibrium 
point q ^ = q 2 = P i = P 2 = 0 anc * takes the form after expansion about this 
point 


H = H„ + H_ + ... +H + ... 
z o m 


12 12 12 , 52 

H 2 = -5 Pi + ? P 2 + q 2 p l ^l p 2 + ‘5 q l ’ kq l q 2“5 q 2 

. . 7JJ k 3 3f3 2 1 1 (3k 2 A 3|I 3 

Ho = - -gj— q, + -ft~ qi ^2 V7 q l q 2 q 2 


q l q 2 +_ TS — q 2 


, 3/ 4 ZDk O 1*0 Z z nr, - ^ n 

*4 = TSf q l + 3T q l q 2 " &T q l q 2 T" q 1 q 2 V!S q 2 


4th k = 


3 {T ( m l ~ m 2 


m^ + m 2 


We now carry out the normalization of the Hamiltonian up to order 3 (N = 3) 
for the Earth-Moon system. For this value of k the eigenvalues corresponding 
to H 2 are distinct and purely imaginary and Assumption 1* holds for N = 3. 
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From (7)^ Hj can be written as 


u _ 1 T - 
Hj - t *■ E r 


where r = (q., q 2 , p ]f p 2 ) and 


E = 


1/4 

-k 

0 

-1 


-k 

- 5/4 

1 

0 


0 

1 

1 

0 


The equations of motion then become 


r = (FE)r + ... 


where 


F = 


0 

0 

-1 

0 


0 

0 

0 

-1 


1 

0 

0 

0 


-1 

0 

0 

1 


0 

1 

0 

0 


As the eigenvalues of FE, X^, X^r -X^, -Xj, are distinct there exists 
an A such that 


A" 1 FEA = D 


( 10 ) 
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where 


D 


0 

0 

0 


0 0 0 

x 2 0 0 

0 -x 1 0 

0 o -x 2 


Moreover A can be chosen so that 


A T FA = F 


(ID 


which guarantees that 

r = AT (12) 

is a eanonieel mapping and thus preserves the Hamiltonian nature of the system. 
As the equations of motion of the system in the new variables become 

r = Dr + ... 


the Hamiltonian in the new variable becomes 

H 2 (q\ / $ 2 ' Pi' ^ = X l Pi + X 2 ^2 ^2 + **• 


Omitting the details, the set of all matrices that diagonalize FE and are 
canonical take the form 


where 




<Vf b !> 


(X 2 a 2 " b 2 ) 


( a 1 ^1 b | ) ( 02 "*” ^bj) 





-1 

6. = A. a, (11-40.), j = 1,2, 

Aj being chosen such that (10) Is satisfied and A so that both (10) and 
(11) are satisfied. The matrix S has the form 



where s^ and Sj are free to vary. As r is real, from (12) it follows that 
we must guarantee that 


Ar = *? 


(13) 
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where the bar represents the complex conjugate. It can be shown that if we choose 
Sj and $2 by 


s 


1 



/ 


S 2 



/ 


then a sufficient condition for (13) to hold is that 






P2 



(14) 


Combining the above matrices, the matrix A has the form 


b/16,1 

(X^-bp/16,1 

(a 1 +A 1 b 1 )/|« 1 | 


1 I 

b 2 / 1 6 2 1 

(A 2 a 2' b 2 )/ 1 6 2 1 

(a 2 +X 2 b 2 V| 6 2 | 


-to/ I 6 1 1 

-i(T^ 1 )/|6 1 l 

-i(a 1 +X 1 b 1 )/|6 1 1 


' °2^ ^ 6 2 ^ 

+i^ 2 /| &2 1 
i(^ 2 a 2 -b 2 )/| 6 2 | 
i(a 2 +X 2 b 2 )/|6 2 l 


and normalizes the terms of second degree of M. 
The Hamiltonian under (12) takes the form 
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~ V l~ v 2~ v 3^ v 4 


H = x i q i p i + x 2 q 2P 2 + 2. VwV 1 ?2 P] P2 


v^+...+ v^ =3 



~ V U V 2_ V 3 *4 


(15) 


Vj+ V2+ Vg+ v^= 4 


v v v v q l q 2 P 1 ^2 + **‘ 
V 2' 3' 4 \ e. \ t 


where (15) is obtained by substituting (12) into (7-9). This relatively simple 
symbolic operation, however, is quite cumbersome when attempted to be done by 
hand. It was done, though, for this particular model with the Earth-Moon con- 
stants and will be used for checking purposes. 


We now normalize the third order terms of (15). As we shall see Assumption 1* 
for N=3 is essential here. Let us introduce a canonical change of variables 
by the contact transformation, [ 3 ] , 


*s» • 

*k " q k + W k 

k = 1,2 

a- _ 3V 
P k \ 


where V (q^, 77 ^) has the form 


V = 


r* ^ v u v 2 v 3 v 4 

° v i ' V 2 #v 3 /V 4 q 1 ^ 7,1 7,2 


v l +v 2 + v 3 + v 4 =3 


(16) 

(17) 


(18) 


We attempt to choose c so as to eliminate as many third-order 

v 1' v 2' v 3 /V 4 

terms as possible in (15). Substituting (16,17) into (15) we obtain 


19 


I> 


- x 1 ^ rj | + x 2 $ 2 n 2 + x 1 - 



+ 



VW V 4 = 3 


r v 2' v 3' v 4 


V 1 v 3 v 4 

| 2 T ij ^2 + terms of degree 4 & higher. 


We note that V is a function of the old variables q£. By a formal process we 

can solve for these variables from (16) and substitute for the q^ in (19). Both 

the transformation from old to new variables and its inverse may be obtained by 

a formal procedure from (16), (17). Both lead to powers series representations 

which converge in a neighborhood of the origin. Eliminating all dependence on 

^ a 

q^ in (19) by this method, H takes the form 


H - X*«+X*«+Xf* *V(e,t?)) 

H - X 1 (■ ] ri ] +\ 2 i 2 ri 2 + \ ] {S ] ‘^1 SrTj > 


(20) 


+ x (* ^(^t?) . iv (4,T?) } + 


v l +v 2 +v 3 +v 4 =3 


g v „ v v in V? nl + terms of degree 4 & higher. 
*1' v 2' v 3' v 4 1 * 1 z 


V 1 v 2 v 3 v 4 

Collecting third order terms in ^ ^2 ^1 ^2 * n and using (18) we 

obtain for a typical term 


c C X . (v. - v«) + X 0 (v- - v .) ] 

v l' v 2' v 3' v 4 3 22 4 


1 /V 2 /V 3 /V 4* 


( 21 ) 


If the bracket in (21) doesn't vanish we can solve for c 

v l'2'3* 

eliminate the corresponding third order term from the Hamiltonian. 


and 

4 

But from 
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Assumption 1* the bracket can vanish only if Vj =v g/ Vj = v^. However 
this would imply the order we are dealing with is even. Thus for N = 3 
Assumption l 1 guarantees' all third order terms of the Hamiltonian can be 
eliminated by the change of variables (16, 17). One must keep in mind 
that although all third order terms are eliminated many more fourth order 
terms arise from this process. These must be kept tract of for error bounds 
and also if higher order normalizations are to be carried out. This has been 
done by hand for the fourth-order terms of the Earth- Moon model and will be 
used for checking purposes. 

Following the same procedure as above, normalization of the Hamiltonian 

up to degree s can be carried out if Assumption l 1 holds for N =s. Let us 

assume that the normalization has been carried out up to degree s - 1 . 

Then the change of variables defined implicitly by (16 - 18) with 

v.+ v„+ v„+ v , = s preserves the normal form up to degree s - 1 . As 
I ^ o 4 Vj Vj Vj v 4 

before, collecting sth order terms in ^2 *1] ^2 leads to an 

equation of the form (21). If jv^-v^i + IVj-v^ j f 0 then 

c can be chosen to eliminate the corresponding sth order term 

V 1, V 2' V 3' V 4 

from the new Hamiltonian. Reasoning as above, all sth order terms can be 
eliminated if s is odd. If s is even then all terms save those for which 
V 1 = v 3' v 2 = v 4' can be eliminated. We choose the corresponding 
c equal to zero in this case. However these terms are formed 

V v 2' v 3' v 4 

from products and lead to an integrable Hamiltonian. It should be noted 

that the complexity of the operation of normalization increases with s. This 
manifests itself in keeping tract of all coefficients that combine to form a partic- 
ular g in (21) which in rum is a function of u!S ptcvious 

V v 2' v 3' v 4 
normalizations. 
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Returning to our original task, after third-order terms have been eliminated, 
the Hamiltonian takes the form 


A 

H ( £]' £ 2 ' ^ \ ' ^2^ = X 1 ^1 ^1 



V 1 + W V 4 



,. v 2' v 3' v 4 1 


+ x 2 i 2 V2 + 

v 2 v 3 v 4 

^2 ^2 + higher order terms. 


Dropping the 4th and higher order terms the differential equations become 


*k X k ^k 

V k = ” X k V 

so that 

«k ■ 

\ = \ (0) e ” X k f . 

It can be shown that if the initial conditions are chosen such that 7)^(0) = • 
then solutions in the original variables will turn out to be real. Thus, if we 
invert all transformation^ information about the original system may be obtained. 
The error in this case comes in because of the truncation of the 4th and higher 
order terms. It is intuitively obvious that this procedure will give better results 
than a linear analysis. For in such a linear analysis the error comes in by 
truncating cubic and higher terms in the Hamiltonian. 

A computer program is being written to perform the algebraic manipulations 
described above. This program is utilizing the I.B.M. FORMAC language and 
is being written for the I.B.M. 7090/94 computer. 
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H. Hermes 


Introduction , The study of "stability* 1 under perturbations, S(t), for a C 
vector field X is no different when the perturbation enters the equation as 
x(t) = X(x(t) + C(t)) , ( x( t ) = or as x(t) = X(x(t))+ C(t) . This is 

no longer true if X is discontinuous. In particular, problems of feedback 
control naturally lead to discontinuous vector fields of the form X(x) = F(x,u(x)) 
where u is a control function. In practice, the value of u is determined 
after making a measurement on the state x(t), at time to If this measure- 
ment is in error, say x(t) + C(t) is measured rather than x(t), the governing 
equation of motion will have the form 

x(t) = X(x(t) + e(t)). (1) 

It is this concept which leads, in section 2, to the definition of stability 
with respect to measurement* Intuitively, X is stable with respect to measure- 
ment if any solutions of eq. (l) and x(t) = X(x(t)), satisfying the same initial 
conditions, remain arbitrarily close over any finite positive time interval 
whenever the supremum of |£(t)| over this time interval is restricted to be 
sufficiently small. 

In general, the initial value problem for a discontinuous vector field 
X need not have a solution. If, however, there is an absolutely continuous 
function cp of the real variable t which satisfies the initial condition and 
cp(t) = X(cp(t)) almost everywhere; we will call cp a classical solution . There 
are many ways to generalize the difinition of solution, so that solutions will 
exist, even if X is merely measurable. A summary of the more standard notions, 
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most of which replace the vector field X by an "averaged" or "smoothed" 
associated vector field, are given by Filippov in [5]. In [3] Filippov defines 
a new concept of a solution, which is tnotivated by control problems; we will 
discuss this notion in §1 and here after refer to such solutions as Filippov 

A 

solutions. 

It will be seen that control laws synthesized from "open loop" controls 
(hence classical solutions exist) may lead to vector fields which are not stable 
with respect to measurement. An example is given for which an optimal feedback 
control exists when solutions are taken in the classical sense, but does not 
exist if solutions are taken in the sense of Filippov. 

The main result shows that if a vector field X is stable with 
respect to measurement (solutions taken in the classical sense in the definition 
of this stability) then every classical solution is a Filippov solution. 

If X is stable with respect to measurement, solutions for t ^ 0 
of the initial value problem for the corresponding differential equation are 
unique, and such a solution when evaluated at a fixed positive time, varies 
continuously with the initial data. This means that, with increasing time, 
solutions may join but not branch. Thus it is felt that feedback controls 
which are meaningful from the viewpoint of applications should lead to vector 
fields which are stable with respect to measurement. To characterize such 
vector fields directly, however, is no easy task. 

$1. A Reason for Fields;, the Filippov Solution. 

Consider a control system of the form 

x = g(x, u(x)) , x = (x 1 ,..,,x n ), U = (u 1# . . . ,u r ) (2) 
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with values u(x) to be chosen from a control set U. Let the terminal manifold 
(target) be a manifold S contained in [0, ») X E n * (E n denotes Euclidean 
n space. ) If g is bounded and Lipschitzian in both arguments and u is a 
given Lipschitzian control, then an initial value problem for (2) with data 
x(0) = x° has a unique solution, with value at time t denoted <p(t,0,x°)o 
Suppose cp(t^, 0, x°) e S. The question considered is the following: If 

S has dimension less than n in E n+ \ is it possible that , for each x in 

some neighborhood yl(x°)<Z 1 E n of x°, there exists a value t(x), 0 ^ t(x) < ®, 

such that cp(t(x), 0, x) € S? 

From a control system viewpoint, it would be desirable that this question 
have an affirmative answer (which is the case if u is allowed discontinuities). 
However for u Lipschitzian we will show, using a method related to a result 
of Br id gland [1, lemma 2], that the answer is negative. Indeed, for fixed t 1 , 
cp(t ! , 0, • ) is a homeomorphism therefore the image of an n neighborhood will 
have dimension Tl. To consider the case where the value of t may depend on 

the point x € /((x°) define the map ^ : E n+1 -> E n+1 by ^(t,x) = ( t ,cp(t ,0,x) ) <> 

Then i|r is a homeomorphism with inverse ^~^(t,x) = ( t ,cp(0,t ,x) ) . Since S 
has dimension less than n, ^~^(S) has dimension less than n. Let P be a 
projection defined by P(t,x) = (0,x). Then P(^"' L (S)) has dimension less 
than no But P(\(r* 1 (S)) is precisely the set of initial points in E* 1 from 
which S is attainable. Indeed x T e P(i _1 (S)) if and only if there exists 
t”£ 0 such that (t f , Cp( t ? , 0, x T ))e S. To see this, x* e P(^ -1 (S))==> for 
some t* , ( t 1 , x* ) e \tr _1 (S) t( t' , x* ) € S or (t 1 , cp( t' , 0, x* ) ) € S. On the 
other hand (t* , cp(t ? , 0, x* ))e S-=>(t f , x f ) = (t ! , cp(o, t* ,cp( t 1 ,0,x* ) ) e ^ ^(S) 
and x f € P(^ _1 (S)). Thus the set of initial points from which S can be 
attained has dimension less than n. 
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It is useful, therefore, from a control viewpoint, to study differential 
equations with discontinuous right sideso For the sake of completeness we will 
briefly discuss the generalized concept of solutions for such equations as given 
by Filippov [3]. 

Let X be a measurable function defined almost everywhere in a domain 
Q E n with values in a bounded set in E n 0 With X(x) associate the convex 
set 

K{X(x)} = n n CO {X(U(x,B) - N)} (5) 

& > 0 (i(N)-O 

where co denotes closed convex hull, U(x, 6) is a closed 6 neighborhood 
of x, N an arbitrary set in E n and p is n dimensional Lebesgue measure. 

An absolutely continuous vector valued function Cp, defined on [0, T], 
is called a. solution in the sense of Filippov of x = X(x) if for almost all 
t, $(t)e K{X(cp(t))}. It is shown in [3] that such solutions will always exist, 
and many of their properties are discussed. In particular, if X is continuous, 
K{X(x) } = X(x) . 

To illustrate this type of solution and its consequences we consider 
a very simple control problem. 

Example 1 . The problem will be that of minimum time transfer with terminal 
manifold S s ((t, x^, x^) : t 1 0, x^ = 0, x^ = 0} and system equations 

*i = 

X 2 = U 2 

with control components subject to the constraint 0 ^ j uj + | u^j ^ 1. It is 
clear that the minimum time needed to attain S from the initial point (x°, x°) 
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is | XjJ + I x^| , and there are many ways in which this can be accomplished. 

We single out two such strategies; each will be given in closed loop (feedback) 
form as symthesized from obvious open loop strategies. 


Strategy 1 



1 

H* 

V* 

O 

y 

x i 

> o , 

O 

All 

x* 

1 

(0, -1) 

y 

x i 

o 

VII 

x 2 > 0 


(1, 0) 

y 

x i 

< 0 , 

x 2 S 0 


(0, 1) 

y 

x i 

o 

All 

x 2 < 0 


l (0, 0) 

y 

x i 

- x 2 

= 0 


Pictorially, the resulting vector field looks as follows: 
Figure 1. 



All vectors are 
unit vectors. 


Strategy 2 


( (o, -1) 


u 2 (x) 


(- 1 , 0 ) 
< <1, o) 
( 0 , 1 ) 


x 2 > 0 

x 2 = 0, x 1 > 0 

x 2 = 


> 


0 , x x > 0 


x~ < 0 



Pictorially: 
Figure 2. 



In each case, the classical solution of the equations of motion exists 
for arbitrary initial data, is uniquely defined for all t ^ 0, depends continuously 
on the initial data and reaches the origin in minimum possible time. These 
same properties are true with strategy 1 when solutions are considered in the 
sense of Filippov, however in the case of strategy 2 the Filippov solutions 
become rest solutions when a state with x 2 = 0 is attained. Therefore, solutions 
in this sense, do net reach the target S. This occurs since the first component 

p 

of the vector field u (x) is zero except on a set of measure zero, i.e. the 
x n axis. From a practical viewpoint, since the control signal is determined 
by a state measurement, one should not expect sets of states having measure 
zero to influence the solution. From this viewpoint, the Filippov solution is 
the more realistic notion. 

In the preceding example, p-^pp-r choice of strategy i.e. 

strategy 1, an optimal feedback control existed whether solutions are taken in 
the sense of Filippov or the classical sense. The following example will show 
that this need not always be the case* i.^. we will produce a feedback control 
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synthesized from optimal open loop controls, which is an optimal feedback 
control if solutions are taken in the classical sense. However an optimal 
feedback control for solutions taken in the sense of Filippov will not exist. 

Example 2. Let the equations of motion be: 


X 1 " x 2 


x p = -x + u , O^u^l , |x°|<2 




with the optimization problem being to minimize the cost functional 
tf 2 2 2 

/ u[(x^-l) + Xg - 1] dt where t^ is the smallest nonnegative time a solution 

o 

reaches the origin. 

2 2 

The open loop strategy of u = 0 until the circle (x^-l) + x^ = 1 

is reached, at which time a switch to u = 1 allowing this circle to be 
traversed in a clockwise fashion, produces a trajectory which reaches the 
origin with zero cost. The corresponding synthesized feedback control leads 
to the following vector field for (4): 


X(x) = 


{ 


On the other hand 


-x 1 + u(x) 
K{X(x 1 ,x 2 )} 


where u(x) 


0 if (x x -l) 2 + x 2 |= 1 

2 2 

1 if (x^-l) + *2 = 1 



since u(x) is 1 only on a set of zero 


measure, and the corresponding Filippov solutions will not reach the origin. 

From the form of the cost functional, it is seen that for any function u(x) 
for which the corresponding solutions in the sense of Filippov reach the origin, 
there will be a positive cost involved. Since this value can be made arbitrarily 
small, but not zero, an optimal feedback control for solutions in the sense of 
Filippov wili not exist. 
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§2* Stability with Respect to Measurement. 


An examination of the two sector fields of example 1 shows a type of 
stability present in the first which is not present in the second. 

For notation ease, if f is a bounded function on a real interval [0,T] 
with values in E n , let ||f|| = ess. sup. {jf(t)| , t € [0, T]}. We will use 
U(x, 6) to denote a compact, spherical neighborhood of radius 5, about the 
point x € E 11 , and coA to denote the convex hull of a set A. 

Definition . A vector field X, for which a classical solution cp of x = X(x) 

with arbitrary initial data x° exists , is said to be stable with respect to 

measurement if given e > 0 and finite T > 0, 3 a & > 0 such that when - 

ever 6 is a measurable function with values in E n and norm less than 6 for 
which a corresponding solution ^ (in classical sense ) of x(t) = X(x(t) + £(t)), 
x(o) — x°, exists on [0, T] , then |Jcp - >|r|j < e. 

For the remainder of this section we will assume X is a measurable 
function defined on a domain Q in E n with values in a hounded set in E n . 

Our concern will be with relating the concepts of stability with respect to 
measurement, Filippov solutions and classical solutions. In particular, lemma 5 
will show that if is a Filippov solution of x = X(x) , x(0) = x° (such 

solutions do exist) then for any e, 6 > 0 , there exists a measurable function 

£ with ||e|| < 5 such that a classical solution cp of x = X(x + 6(t)) , x(0) = x° 
exists and satisfies ||q> - >|r|| < e. This essentially says that if one allows 
arbitrarily small perturbations of the argument, a response to any vector field 

X may he made to agree closely with a response to the associated Filippov 

generalized field K{X(*)}. After this has been established one easily obtains: 

Theorem 1. If X is stable with respect to measurement then every classical 
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solution is a Filippov solution. 


Lemma 1 . Let if be an absolutely continuous function on [0, T] vith values 

in E n , and z a measurable function with z(t)e K{X(t( t) ) } , t € [0, T]. Then 

1 k 

given any 6 > 0, there exist a finite number of measurable functions u) , . • . ,u) 
with ab^tje U(*(t), 6) such that the function v 1 defined by v 1 (t) = X(o) 1 (t)) 
are measurable and for any e > 0, z(t) i^s contained in an e neighborhood of 
co(v 1 (t) , . . . ,v k (t)) . 

Proof Filippov [3] shows that the requirement z(t)e K{X(\|r(t))} is equivalent 
to the condition that for any vector t] , 

z(t)-T) ^ lim (ess. max (X(u)‘T) : u€U(i|r(t), r))) (5) 

X -»0 

or equivalently z(t)-T) ^ ess. max {X(u)*r] : ueU(^(t), r)} for every f > 0. 

Let z be any measurable function with z(t)e K(X(\|r(t ) ) } . Suppose we 
are given 5, € > 0. Pick an arbitrary vector tj 0 ; we will first show that 
one can construct a measurable function cr with o>(t)e U(i(t), 5) such that 

z(t)'T) § X(cD(t))«rj + e/2 (6) 


for t e [0, T]. 

Subdivide the interval [0, T] into subintervals by a partition 

0 = t < t_ < . . . < t = T and let 5(t-t.) be a continuous real valued function 
o 1 m i' 

defined on [ t^, t^ ) with 6(0) = 5, 6 ^ 5(t-t^) ^ 5/2 and such that 

u(i(t), 6(t-t.)) c u(+(t’), s(t'-t.)) (7) 

for t. St' S t < t, , . The existence of such a partition and function &(t-t. ) 
i l+l i 

is an immediate consequence of the uniform continunity of i|f on [0, T]. 
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Figure 4 



For any > 0 , by Luzins' theorem we may cloose a compact subset 

Ei(e ) of U(t(t^), &) differing in measure from U(+(t^), 6) by less than 

and such that X is continuous on E^(c^). Then for every t e [t^, t^ + ^) 

the set = E^(e^) ^ U(i(t), 6(t-t^)) is a nonempty compact subset of E^(e^) 

on which X(*)»rj is continuous. Let v (t) = max.{X(u>)*i} : co € K + ). By (7), 

€ 1 1 

v € is a monotone decreasing function on [tj, t^ +1 ) hence measurable. By 

theorem 1, [2] there exists a measurable function co with u>(t) € such that 

X(cu(t) ) • T) = v e (t) , t e [t^, t^ + ^). (Here we have replaced the condition of 

the sets expanding, C- K t , for t < t' in the cited theorem, by 

contracting, K , C- K , for t' < t ; a condition which does not alter the proof 
x> t 

since the direction of traversing the time axis is immaterial.) 

This defines the function cd on the subinterval [t^, t^ + ^) j since 
i was arbitrary we may assume cn to be defined on [0, T] as that function 
whose restriction to [t.,, t., Al ) is defined as above. 

For any €.. > 0 either v (t) S ess. max. (X(u)*T) : u € U(t(t), B/2)) S 

1 € 1 

z(t)»T| > the latter inequality following from (5) , or v (t) < ess max. {X(u) *tj : 

1 

u e U(t(t), b/2)}. In the first case we deal with the situation where the 
maximum of X(*)'tj occurs in U(t(t), 6)- U(+(t), 6/2) and inequality (6) holds 
even with e = 0. 
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In the second case, for fixed t, v (t) is an increasing function 

e l 

°f g^ since we may assume, without loss of generality, that E^(g^) 

if e^ < e^. Since, in this case v^ (t) < ess.max{X( u) • r\ • u e U(f( t) , 5/2) } , 

we may take a sequence of values e^ tending to zero. The corresponding 

bounded monotone sequence of real numbers v (t) must converge to ess. max{ X( u) • t) 

€ i 

u € U(^(t), 5/2)} ; indeed if it converges to a number m less than this, by 
the definition of ess. max. there exists a set F of positive measure on which 
X( • ) • T] > in . To obtain a contradiction we need only choose e^ less than the 
measure of F. 

This establishes that for e^ sufficiently small, there exists a 
measurable function oo with values cjo(t)e U(i(t), 6) such that v(t) = X(<n(t)) 
is measurable and inequality (6) is satisfied. (Note; It is not true in general 
that a measurable function of a measurable function is measurable.) 

Now let S n ~^ be an n-1 sphere in E n which contains 

^ n 1 

U +€ [ 0 T jX(U(\lf(t), 6); this exists by hypothesis. Since S is compact 

choose a finite number of vectors T) 1 , i = 1, 2, ...,k belonging to S 

and so that e /2 neighborhoods of the rj 1 cover S n_1 . For each T] 1 construct, 

as before, a function go 1 , measurable with values of^t) GX(U(^(t), 5)) satisfying 

(6). Let v 1 be the corresponding measurable function; v^(t) = X(o) 1 (t)). Then 

1 k 

z(t) is contained in an € neighborhood of the convex hull of (v (t),...,v (t)}. 

1 k 

Lemma 2 . Let v , . . . , v be bounded measurable functions defined on [ 0, T] with 
c 4(t) = { v X (t),. . .,v k (t)) . Let co v/ 4( t) denote the convex hull of ^(t). Then 
if z is a measurable function with values z(t) contained in an € neighborhood 

t 

of co^(t), there exists a measurable function v with values in c/4(t) such 
t ~~ 

that |/[z(t) - v(t)]cK| < e(T+l) uniformly for t g [0, T]. 
o 

Proof . a) We will first show there is a measurable function y with values 
y(t)e cor^(t) such that ||z-yj| ^ g. 
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Using the terminology of m, if z(-) is measurable single valued 

function then U(z( - )> c) is a measurable many valued function. Indeed, if 

B(y°, r) is a closed ball of radius r, center y°, {t : U(z(t), e) fl B(y°, r) (j)} 

= { t : j z(t) - y°| £ r + e} which is measurable. 

Next, since the functions v 1 are measurable, we will show coc^(*) 

is a measurable set valued function. Obviously co c >^(t) is nonempty and closed 

o n 1 

for each t. Letting B(y , r) be as above and S denote the unit n-1 

sphere we note that the distance from co{ ( v^(t )-y°) , — , ( v^(t)-y°) } to the 

origin Is max (min tj» ( v X ( t)-y°) ) . Then 
T]eS n-1 ISiSk 

{t : co cA{t) D B(y°, r) =f: (J)} = {t : max (min tj* (v^(t)-y°)) £ r} 

tjcS 11 " 1 lSi^k 

which is measurable. 

From [4], U(z(*)> c) fl co • ) is again a measurable set valued 

function and there exists a measurable single valued function y with 
y(t)e U(z(t), e) fl coo^(t). 

b) We next show that if y is a measurable function on [0, t 1 ] 

with y(t)e co o4{t) for each t € [0, t 1 ] then y admits the representation 

y(t) = Z a i (t)v 1 (t) where the scalar valued functions are measurable, 

k 

0 £ a^Ct) £ 1 and Z^^a^t) = 1 for all t € [0, t 1 ]. 

This result is closely related to lemma 1 [6]; which would yield the 
desired result if the functions v 1 were continuous. To modify this to the 

i k i 

present case where the v are measurable, let f (t, a) ~ or/ 1 ' (+) ; 

Q = {a € E k t =1, 0 S a i § 1}, and R(t) = f (t, Q,). Then f if 

continuous in a for each fixed t. Referring now to the proof of the lemma 
of Filippov [5] and letting QL, v 1 play the role of the u^, z* , respectively. 
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1 k 

in that, proof, choose the set E to he so that ...,a , ...,v and 

y are continuous in E. Because of the special form of f it follows that 
f(t, a) is continuous on E X Q and the Filippov argument may be applied 
to give the desired representation for y. 

c) From theorem 1 [6], it now follows that for any interval [0, t 1 ] 

there exists a measurable function v with values v(t)e^>^(t) for each 

t € [0, t' ] such that 

t* t* 

/ y(x)dx = / v(x)dT . (8) 

o o 

Since the functions v 1 were bounded there is a constant M such 
that ||y|| ^ M , || v|| i M. Subdivide the interval [0, T] into m equal sub- 
intervals each of length T/m. Let I. denote the interval (jT/m, (j+l)T/m]. 

Using (8) for each j = 0, ...,(m-l) , define v on I . so that / T [ u( t) -v( t) ] dx = 

J J t 

0. Now if m is chosen so large that m i 2MT/e it follows that | / [ y( t) -v( t) ]dx| 
< e uniformly for t e [0, T]. 

d) To finish the proof we show the function v constructed in part 

t 

c) satisfies the conclusions of lemma 2, Indeed |/ [z(t)-v(T)]dT | = 
t tot 

I / [z( T )-y(T) + y(T)-v(T)]dT I S 1/ [ z( T)-y( t) ]dr| + | / [ y( t)-v( t) ]dx| £ 
o o 

cT+€=€[T+l]j using the results of a) and c) , respectively. 

Lemma 3 . Let be a Filippov solution of x = X(x) , x(0) = x°. Then for 
any e , 6 > 0 there exists a measurable function £ : [0, T] -» E n with 
||c|| < B such that a classical solution cp exists , on the interval [0, T], 
for the problem x = X(x + e(t)) , x(0) = x° and satisfies ||<p-\|r|| < e. 

Proof . Let i)f( t) = z(t)e K{X(\|f(t) )) . By lemma 1 there exist k measurable 

1 k i i 

functions <n with ad (t)e U(>|r(t), 6/2) such that the functions v 
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defined by v 1 (t) = X(o) 1 (t)) are measurable and z(t) is contained in an 

e/(T+l) neighborhood of co{ v^ft), . . . , v k (t) } . 

By lemma 2, there exists a measurable function v with v(t)e(v^(t), 
k t 

. ..,v (t)} such that j/ [v(x) - x(x)]dx| < € uniformly for t e [0, T]. 

o 

Define the measurable function cu by Go(t) = go (t) if t is such that 

v(t)=v i (t). Then go is measurable, oo(t) e U(^(t), 6/2) and X(o>(t)) = v(t). 

We next produce the absolutely continuous function cp and measurable 

function 8 in the statement of the lemma. 

t t 

Define qp(t) = x° + / v( r)dn Then | c P(t)-t(t)| = j/ [ v(x)-z( x)]dxj < e 

o o 

for t e [0, T] hence cp(t) e U(t(t), e) and ||cp-i|r|| ^ e. Define 8(t) = 

go( t ) - cp ( t ) ♦ Then 8 is certainly measurable and |e(t)| = | oo(t) -♦(t)+>lr(t)-cp(t) | ^ 

^ 6/2 + e. There is no loss in generality if it is assumed € < 5/2. Therefore 

II e|| < 6. 

Also, cp(t) + S(t) = co( t ) hence X(cp(t) + 8(t)) = v(t) and from the 

t 

definition of tp, cp(t) = x° + / X(cp(x) + 8(x))dx showing that <p is a classical 

o 

solution of x = X(x + 8(t)), x(0) = x°. 

Proof of Theorem 1 : We shall prove the contrapositive j i.e. if some classical 

solution exists and is not a Filippov solution then the field X is not stable 
with respect to measurement. 

The assumption that some classical solution is not a Filippov solution 

o o 

implies there exists x and classical solution qp through x existing on 

some interval [0, t n ] such that there is a Filippov solution f through x° 

with cp(T) - ♦(T) 4 0 for some T e (0, t^] . Let |cp(T) - y(T)| = r > 0, pick 

e = r/2. Then by lemma 3, for |je|| arbitrarily small, we can find a classical 

solution | of x = X(x + 8(t)), x(0) = x° such that | £(T) - +(T)| < e 

hence | cp(T) - £(T)| > e, i.e. X is not stable with respect to measurement. 
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We shall end by briefly summarizing some additional properties of 
vector fields which are stable with respect to measurement. 

If X is stable with respect to measurement, the solutions of 
x = X(x) for t ^ 0, are unique. This follows as an immediate consequence 
of the definition. 

If X is stable with respect to measurement and cp(t^> x°) denotes 
the solution through initial data x° evaluated at time t-^ > 0, then cp(t^,*) 
is continuous. 

k. o 1c o 

Indeed suppose x x but cp(t^, x ) cp(t^, x ). Then there 

exists a 6 > 0 such that Icp^, x k ) - cp(t , x°)| §5 for all k sufficiently 
k ok 

large. Let 6 (t) = x - x ; i.e. a constant measurement error. For k 
sufficiently large, ||e|| can be made arbitrarily small. 

Since <j>(t, x k ) = X(cp(t, x k )) , if we define | k (t) = cp(t,x k )-x° + x k 
then | k (t) = <p(t, x k ) hence £ k (t) = X(| k (t) + G k (t)) and | k (0) = x°. From 
the definition of £ k , for k sufficiently large || | k - cp(’, x k )|| can be 
made arbitrarily small j it follows that || (; k - cp(*, x°)|| > &/2 for k suf- 
ficiently large, hence X is not stable with respect to measurement. 
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ABSTRACT 


A survey of various approaches to the problem of estimating 
the domain of attraction of an equilibrium solution to a system 
of nonlinear autonomous differential equations is given. Based 
upon observations resulting from this survey the problem is re- 
formulated as that of optimally choosing the Liapunov function 
from the space of positive definite quadratic forms . An esti- 
mate of the domain of attraction is then obtained as the solution 
of a minimization problem. This approach to the problem has 
the advantages of: 1) being designed specifically for machine 

computation; 2) yielding an estimate that is readily visualized; 
and 3) being relatively insensitive to system dimension. Some 
preliminary numerical results are presented for the Duffing 
equation with damping. 
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Introduction 


This presentation is concerned with the problem of computing 
the restrictions on the initial state errors in a dynamic system 
which will guarantee that these state errors will tend to zero 
as t-*». Thus, we shall be concerned with developing an effi- 
cient numerical technique for estimating the domain of asymptotic 
stability or more succinctly the domain of attraction of the 
null solution. This problem has application in qualitatively 
predicting the attitude motions of a satellite and perhaps may 
provide a first step toward solving the problem of qualitatively 
evaluating the effects of disturbances upon various rocket 
guidance schemes . 

The applicability of this analysis to rocket guidance prob- 
lems is crucially dependant upon the assumption that motions of 
the vehicle off the nominal path can be described by an autonomous 
state differential equation, viz., 

x = g(x,u(x>) = h(x), h(0) = 0 (1) 

where x(t) is the n- vector describing the deviation from the 
nominal state, u(x) represents the control lav? designed to 
control this deviation, and the null solution is an equilibrium 
solution. The domain of attraction Q is then defined as the 
set of all initial points that generate trajectories that tend 
toward the equilibrium solution, i.e., 

a : (x° | ^im x(t;x°,t Q ) = o) (2) 


The only body of theory that has been applied to the general 


_ £Z 
OJL 
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pro Diem or estimating uic uOiuaUi 
direct method. Within this theory there are two distinct ap- 
proaches that have been taken to determine the domain of attraction. 


The first of these approaches, due to V. I. Zubov [1], 
allows an exact solution to the problem, if an arbitrary function 
can be chosen such that a closed form solution is obtained lor 
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Zubov's partial differential equation, or it allows an estimate 
of the domain of attraction via a truncated power series solu- 
tion to the equation. That is, if a positive definite 0(x) 
can be found such that the Zubov linear partial differential 
equations 


V x v(x) • h(x) = -0(x) (l - v (x)) 
or 

V x v(x) * h(x) = -6(x) (l - v (x)) (l + h(x) » h(x)) 
can be solved exactly for v(x) , then Cl is given by 

Cl : (x | 0 < v(x) < l) 

If a power series solution is obtained in the form 

n ^ 

v n (x) = £ v i( x ) » v i( ax ) “ a v t (x) , 
i=2 

a series of homogeneous forms, then an estimate Cl of Cl 
is obtained via 


ci 


n 



0 < v n (x) < l) 


> 


and 

ci n C ci . 

In 1962 Margolis and Vogt [2] reported on a procedure which 
employs a digital computer to develop the series solution to 
Zubov's equation for a class of differential equations of dimen- 
sion two. The authors noted two principal problems: 1) com- 

putational problems arise for systems of higher dimension; and 
2) the convergence of the series solution is far from uniform, 
i.e., the estimate obtained for the Van derPol equation by 



using only second degree terms in the series was better than the 
estimate obtained by including terms up to sixth degree. 

During 1962 and 1963 Szego [3] reported on methods for 
solving Zubov's equation in vector-matrix form (these methods 
are related to his earlier work [4] on generating Liapunov 
functions via a "quadratic form" whose coefficients are functions 
of the state variables) and in [5] on a generalization of Zubov's 
equations . The latter was pursued somewhat further by Szego and 
Geiss [6]. Although some results regarding identification of 
limit cycles and estimation of domains of attraction are reported, 
no results regarding the conversion of these processes to numer- 
ical algorithms are given. 

Rodden [7] and [8] reported in 1964 on an algorithm he 
developed for both calculating the solution to Zubov's equation 
in series form and analyzing the resulting Liapunov function. 

His work was restricted to problems of dimension two and three, 
and his results indicated three principal problems: 1) lack of 

uniform convergence of the series solution to Zubov's equation; 

2) strong dependence of the final result upon the choice of the 
arbitrary or "constraint" function 0(x) in Zubov's equations; 
and 3) visualization of the estimate of the domain of attraction, 
particularly for three dimensional systems. Rodden found, in 
some examples, that the second degree approximation was better 
than the 20th, and that the convergence of this series solution 
could be improved by solving a modified Zubov equation. How- 
ever, this change still did not guarantee that higher order 
approximations would be better than lower order approximations . 

The second principal approach to estimating the domain of 
attraction is to base the analysis upon La Salle's theorems on 
the extent of asymptotic stability [9] and use one of the many 
procedures for developing Liapunov functions that are available 
in the literature [10] and [11]. This tack was reported on in 
1962 by Infante [12] and Infante and Clark [13]. Infante developed 
an ingenious and successful procedure for developing Liapunov 
functions for two dimensional systems based upon an approximation 
to the dynamic system. Although estimates are easily obtained 
from his Liapunov functions, the technique for generating the 
functions does not appear to be suited to machine computation. 
Infante's work was developed in 1964 by Walker [14] for systems 
of higher dimension but again the technique is not suited to 
machine computation. The present author [15] reported in 1964 
some favorable results obtained from a cursory look at the value 
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of using "optimal" quadratic form Liapunov functions for esti- 
mating the domain of attraction. (This paper reports on an 
extension of this concept.) In 1965, Weissenberger [16] and [17], 
using the analysis algorithm developed by Rodden [7], developed 
a numerical technique for estimating the domain of attraction 
of relay control systems via an "optimal" choice from the class 
of Lure -Liapunov functions . 

Thus , upon reviewing the history of this problem the follow- 
ing remarks become apparent: 

1. The majority of techniques for generating Liapunov 
functions are unsuitable as bases for machine 
computation of Liapunov functions because of 

the requirement of experience and ingenuity in 
their application. Of those which are acceptable, 
i.e., series solution of Zubov equations, Lure- 
Liapunov formulation, and quadratic forms, the 
Zubov approach suffers from erratic convergence 
and lack of knowledge of how to choose the 
"constraint" function. 

2. The method of analyzing the Liapunov function 
to determine an estimate of the domain of 
attraction should be relatively insensitive 

to system dimension. Rodden 's technique depends 
on geometric analysis to determine points of 
tangency of hypersurfaces and hence is directly 
dependent on system dimension. 

3. The estimate of the domain should be easy to 
visualize if it is to be of engineering value. 

A glance at the figures constructed by Rodden 
for three dimensional problems, and recognition 
of the fact that rocket guidance systems are 

of at least dimension four gives strong motivation 
to this statement. 

4. Little attention has been given to selecting 
the "optimal" Liapunov function from a given 
class of functions; rather, the emphasis has 
been on new methods of generating Liapunov 
functions . 
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Based upon these remarks, the Liapunov function to be used in 
this analysis will be restricted to be a member of the class of 
positive definite quadratic forms. This restriction guarantees 
that the estimate of the domain of attraction will always be an 
ellipsoid and thus easier to visualize than the results of higher 
order estimates. Secondly, based upon the results of Margolis 
and Vogt [2], and Rodden [7], there is reason to believe that 
this estimate may be better than those obtained by using functions 
of higher degree, particulary if the quadratic form parameters 
are optimally chosen. Finally, information may be gained that 
will aid in formulating a best choice of the "constraint" function 
0(x) for Zubov's equations. 

Problem Formulation 


Consider the basis of this analysis, i.e., 

Theorem (La Salle [15]): 

Let V(x) be a scalar function with continuous first 
partial derivatives. Let £2^ designate the region where 
V(x) < l. Assume that £2^ is bounded and that within £2^: 

V(x) > 0 for x M 

V(x) < 0 for x £ 0 , 

Then the origin is asymptotically stable, and above all, 
every solution in £2^ tends to the origin as t-*-oo. 

Thus, £2^ is an estimate of £2 and the problem is reduced 
to choosing V(x) from the class of quadratic forms and then 
establishing that the required properties exist in some domain. 
That is, the following must be accomplished: 

1. Prove that V(x) is positive in some region 
that includes the origin. 

2. Prove that V(x) = W • h(x) is negative in 
some region including the origin. 

3. Establish a region £2, within which both 1 and 

2 hold . 1 

4. Prove that £2, is bounded. 

I 
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Now since V(x) is restricted to be a positive definite quadratic 
form, viz., 


V(x) = x T Px , P > 0 (9) 


it is positive everywhere. Further, restrict the system equation 
(1) to have a stable linear part, i.e., let 

x = h(x) = Ax + f(x) (10) 


where A is a stable matrix (its eigenvalues all have negative 
real parts) and f(x) contains no terms of first order in x. 

This is not an inordinate assumption since our present technology 
only allows synthesis based upon essentially linear analysis and 
thus a system with stable linear approximation usually results . 

# Based upon the assumption of equation (10) the calculation 
of V(x) results in 

V(x) = x T (A T P + PA)x + 2x T Pf(x) (11) 


and choosing 

- Q = A T P + PA (12) 


results in 

V(x) = -x T Qx + 2x T Pf(x) . (13) 

Thus, if Q is positive definite V(x) will be negative in a 
region including the origin by virtue of the fact that f(x) 
contains no terms of first order in x. Now, since A is 
assumed stable we know that every positive definite Q will 
produce a positive definite P (Kalman and Bertram [18]) and 
thus requirements 1 and 2 are satisfied. 
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The establishment of Cl. is next on the agenda. By virtue 
of our restriction on V(x) the set of hypersurfaces V = c£, 

C 1 < c 2 < c 3 < • • • < c i < • • • < cn will be a set of ellip- 
soids of fixed orientation and increasing size. Thus Cl ^ 
should be chosen to be the interior of the largest such ellip- 
soid within which V < 0, and that ellipsoid will be the 
smallest one which has a point of contact with the hypersurface 
given by $ = 0, x ^ 0 (see Fig. 1) . 



Fig. 1 Typical Relationship of the Loci V = 0 and V = constant 


53 





The problem is then one of calculating 


f, = min V(x) subject to V(x) =0 . (14) 

x -^0 

It is exactly at this point that a computer is of most value . 

Also, note that since 

: (x | V(x) < i'j (15) 

is an ellipsoid it is bounded and requirement 4 is satisfied. 

We have tacitly assumed above that there is a solution 
to the problem stated in (14) , a sufficient condition for exis- 
tence is obtained as follows . Consider that we must prove that 

-x^Qx + 2x^P f(x) < 0 (16) 


in some domain D including the origin. Now note that 

x T Pf(x) < |x T Pf(x)| £ ||x T P|| | | f (x) | | (17) 


via the Schwartz inequality, and using the extremal properties 
of characteristic values of pencils of quadratic forms, 
Gantmacher [ 19 ] , 


||x T r||< vA-V) INI - *““<p> 1UII 


(18) 


where A (P ) is the maximal eigenvalue of P . Similarly, 


x T qx > ? mln (q) |UI| 2 • 


( 19 ) 


,min 


where A (Q) is the minimal eigenvalue of Q. Thus, (16) 
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is satisfied if 


i|f(x)|| < 


2^ mn (P) 


II: 


in 


( 20 ) 


Since Q and P are positive definite 


,min 

A 

2A max 


£21 

(P) 



> o 


( 21 ) 


and thus (20) is a special case of a Lipschitz condition. 

The results of this procedure are dependent upon the choice 
of Q and for each Q there will be a different , thus 

perhaps there is a best choice of Q for a particular criterion 
function. The most obvious criterion is the volume of 
and thus the last step in the analysis is to define 


j ( q ) - ^ /2 ( n \( p ))' i/2 

i=i x 


( 22 ) 


and 



J(Q°) = max J(Q) 
Q>0 


The 

computational procedure will then 

be as follows: 

1. 

choose Q > 0 


2. 

calculate P via (12) 


3. 

compute £ via (14) 


4. 

calculate J(Q) via (22) 


5. 

modify Q in direction of larger 

J(Q) 

6. 

return to 2 



(23) 
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Numerical Results 


Consider the Duffing equation with damping, viz.. 


X 1 = x 2 

• , 3 

x 2 = “ e l x 2 “ X 1 + e 2 X l 


(24) 


and the quadratic form Liapunov function 


V = p u xj + 2p 12Xl x 2 + P 22 * 2 


(25) 


whose time derivative with respect to (24) is 


V = - 


2 2 
2p 12 X l ^ 2p 22* 1^12 "* X 1 X 2 ^ 2e l P 22 “ ^P]_2^^2 


2e 2 p 22 x l x 2 + 2s 2 P 12 x l 


( 26 ) 


Thus, the following relationships exist: 


A = 


P = 


Q = 


0 1 \ 

, f(x) = 

1 1 

(°, ) 
\ e 2 X l / 

P 11 P 12 \ 


P 12 P 22 ' 


2 p 12 

P 22 + e l P 12 

P 22 + € 1 P 12 " P 11 

2e l P 22 


12 


(27) 
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Application of the Sylvester criterion for positive definiteness 
of Q and P yields the parameter restrictions (for = 1) : 


( 


P 12 




P 12 > o 


( 28 ) 




Pll > 0 » P 22 > 0 


(29) 


This system has three equilibrium points , viz . , 


(( 0 , 0 ) 


(x 1 ,x 2 ) 


<(^7»o) 

v (V^2*°) 


(30) 


the first being stable and the others unstable. Thus, one would 
not expect the domain of attraction to exceed |x^| = £2 ' 2 
and it is reasonable to inquire whether (20) is satisfied in 
D, where 


D : (x| | |x| | 2 < e 2 X ) 


Hence, the question is what is 


K 


such that 


(31) 


= 2^1 , w < K 2 in D 
1 |X1 1 V^i+ x 2 


( 32 ) 


and the solution is 
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1 


(33) 


e 2 X l 


K > SU P 7 T=l. 

v /x 1 + X 2 


xeD 


Thus, the estimate of the domain of attraction will be larger 
than the circle ||x|p =e"^ if 


K 2 = 


A mln (Q) 

2A max (P) 


> 1 


(34) 


or if 


p -/7 p~ 


2) + (P - a + 1) 


> 4 


(35) 


where a = Pn(Pi 2 ) _ ^ and ^ = P22(P12) _i - The parameter 
restrictions (28), (29) and (35) are illustrated in Fig. 2. 

The allowable choice of parameters is as given in (28) , 
(29) and (35) or Fig. 2 and the analysis in [15] has shown 
that the optimal choice, using the area of as criterion, 

to be 1 


,-l 


'll 

5 12 


= 2 


122 

>12 


= 1 


which is exactly on the boundary of the allowable region. Note 
that the resulting V is semidefinite, i.e., V = 0 on = 0, 
and one must use another form of the stated theorem. See [15], 
or [9], p. 66. The corresponding (for e£ = 0.04) is 

shown in Fig. 3 along with the estimate obtained by using the 
energy of the undamped system as the Liapunov function. The 


(36) 
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Fig. 3 


Comparison of Best Estimate, 0^ , with Energy Function 
Estimate, 0 . , and System Trajectories 
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value of i is 12.5 when p-^ = 4 and thus the estimate is 


n o 1 2 1 

n t : 2 X 1 + 2 x l x 2 


1 2 

+ £ < 12.5 


( 37 ) 


The points of contact with ^ = 0 are at the equilibrium points 
(xi, X 2 ) = (+5, 0) and the estimate is seen to be close to 
the actual seperatrix and considerably larger than the estimate 


: \ K i + \ x 2 “ 0«01x^ < 6.25 
< 25 

obtained by using the energy of the undamped system as the 
Liapunov function. 


(38) 


In Fig. 4 the optimal quadratic form estimate, fig , is 
compared with the estimates obtained by Infante [12] for the 
system with €2 = 1. Ingwerson's procedure [20] for generating 
Liapunov functions yields the following estimates: 


4 2 

2 X 1 x 2 1 

fl l : X 1 4 + X 1 X 2 + 2 < 4 


(39) 


while Infante's procedure yields 


: x i 


2 
X 1 

2 

Z 2 


2 3 

+ 2xjX 2 + x 2 < j 


lxj < 1 


(40) 


Thus for this example it seems that the optimal quadratic form 
technique yields an improvement in accuracy, ease of solution, 
and ease of portrayal of the estimate of the domain of attraction. 


The numerical solution of the constrained minimum problem 
(14) is obtained by solving an unconstrained problem, viz.. 
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Fig. 4 Comparison of Best Estimate, , with Estimates 
Obtained by Infante, and System Trajectories 


6 2 


V(x) + 


(41) 


l = min 
x 



* 


2 

where 2 ll x ll was introduced to avoid the trivial solution 
and K| is manipulated to assure satisfaction, to a prescribed 
accuracy, of the constraint V(x) = 0. 

At present, we are using a new algorithm for finding the 
minimum of a very general class of functionals to solve (41) . 
(Eventually it will also be used to solve [23] .) This 
algorithm is being developed at Grumman by Mr. R. McGill based 
on work by Davidon [21]. The algorithm being developed by 
McGill has the following salient characteristics: 

1. It does not require numerical inversion of linear 
operators and thus is relatively free of 
dimensional limitations. 

2. It is stable with respect to convergence, i.e., 
convergence to a local minimum is guaranteed. 

3. It is efficient, i.e., convergence is quadratic 
in a neighborhood of a minimum. 

4. It allows a tradeoff between precision and 
computing time . 

5. It requires modest storage. 

Typical computational results for G 1 = 1, e 2 = 004 
are presented in Figs. 5 through 12. These figures show the 
boundary of the estimate and its relationship to the 

constraint V = 0, and the area contained with P ^ . The esti- 
mate Q.£ is the elliptical region surrounding the origin. 

The aberrations from ellipticity are due to the plotting machine 
routine being used emu are not part of P.£ . Th^ inn' v = 0 
are those with the triangle and square markings . These mark- 
ings are used to distinguish the branches corresponding to the 
positive and negative roots, of the quadratic equation (in X£ ) 
used to generate the loci V = 0. These markings along the 
x^ axis indicate that the roots are complex for those values 
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of and are not part of the loci V = 0. Note that when 

PH = 4.0, pi2 = 0.5, P 22 = 2.0 V = 0 appears to have a cusp 
at its point of contact with V = i . (This situation would be 
difficult to handle using a geometric approach such as Rodden's). 
The best computed result is about 107 o off the optimal 
j(qO) = 50-ir = 157. Some convergence difficulties have been 
observed as the boundary 2 of Fig. 2 is approached. This 
phenomenon has not yet been investigated. 

Conclusions 


Results obtained by investigators who pursued estimation 
of the domain of attraction via Zubov's technique, and the 
preliminary results presented here indicate that estimation 
of the domain of attraction for quasi linear dynamic systems 
via "optimal" quadratic form Liapunov functions, as formulated 
here, is feasible. This procedure, in conjunction with McGill's 
algorithm, offers the advantages of: 1) readily leading to an 

efficient algorithm for estimating the domain of attraction 
which is relatively insensitive to system dimension; and 2) 
providing an estimate which is easy to visualize, i.e., an 
ellipsoid. 
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(7 


In the study of trajectory optimization by classical calculus of 
variations techniques, one normally encounters certain Euler-La grange 
equations involving the control angles used in the formulation of the 
problem* In some cases, these equations lead to a solution for the 
angles (or, usually, the tangents of the angles) in terms of the 
Lagrange multipliers; and these solutions can be used to eliminate the 
control angles from the Euler -La grange equations and thus leave a sys- 
tem of differential equations in the state variables and the Lagrange 
multipliers* This is generally desirable. However, it seems that the 
process is more readily carried out in some coordinate systems than 
others and, in fact, virtually impossible in some systems. Thus, if 
one happens to be using a coordinate system in which the latter is true, 
it would be convenient to transform to another coordinate system in 
which the problem was not present, find the desired solutions, and then 
transform back to the original system* This involves transforming the 
state variables and their corresponding Lagrange multipliers from one 
system to another. In this discussion, the technique for a general 
transformation of this type is given and then applied to a specific 

problem involving three-dimensional trajectory optimization in a plumb- 

• * • 

line coordinate system (with state variables x, y, z, x, y, z, and 
control angles X and X and in a spherical coordinate 

pitch y aw 

1 

W. E. Miner, "Methods for Trajectory Computation", NASA-MSFC, 
Aeroballistics Internal Note No. 3-61, May 10, 1961 
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system (with state variables r, <t> , 0 , v t Y , S an d control angles 
2 

a and 8 ) , In the former system, tan X . and tan X are 

pitch yaw 

readily solved for whereas the same is not true in the latter system* 
However, the desired result is obtained for the latter system by appli- 
cation of the method just outlined. 

Consider the equations of motion which simulate vehicle flight 
in three dimensions, through a vacuum, for a non-rotating spherical 
reference body. Thrust and weight flow are assumed constant and thrust 
and gravity are the only two forces acting on the vehicle. These 
equations, in flight path coordinates, are: 


2 

D. H. Young, "Three Dimensional Vacuum Trajectory Optimization 
with End Points in Flight Path Coordinates", Douglas Aircraft 
Company Space and Missile Systems Division Memorandum 

Symbols from the above are as follows: 

v: total missile velocity, directed along the flight path 

Y: vehicle elevation flight path angle (angle between the 

projection of the velocity vector on the local tangent 
plane and the velocity vector) 

5 : vehicle azimuth flight path angle (angle between north 

and the projection of the velocity vector on the local 
tangent plane- positive, clockwise from north) 

a : in-plane angle of attack (angle between the velocity 

vector and the projection of the thrust vector in the 
v-n plane- positive, counterclockwise from the velocity 
vector ) 

8 : out-of-plane angle of attack (angle between the velocity 

vcuLui and the projection thp thrust vector in the 
v-s plane- positive, clockwise from the velocity vector) 
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V 


- g sin Y 


T tan a 


Y = ■ . . . — . - + — 1 cos Y 

Mv( i Vtan^ a + tan 6 + 1) L r v 


T tan 8 


Mv cos y ( Vtan 2 a + tan 3 + l7 


+ v cos Y sin 5 


r cot 


r = v sin Y 


<t> = (v cos Y cos 5 )/r 


0 = ( — v cos Y sin <$ )/r cos <P 


where r, ji, 9 are the usual spherical coordinates and v, Y , <$ f 
ot , and 8 are as previously defined. 

In order to use classical calculus of variations techniques, 

form: 


i 

.M( ± V^tan^ a + tan^ 3 + l) 


- g sin Y 


XY 


T tan a + (l _ l\ 

Mv( * /tan 2 a + tan^ 8 + 1) \ r v / 


cos Y 


1/ X £ 


T tan 


L Mv cos y ( 1 \/tan^ a + tan 2 8 + 


+ 


v cos Y sin 


1) r cot 4> 
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A [v sin y ] + X* [(v cos y cos 6 )/r] + 

r <P 

X Q C(-v cos y sin 6 )/r cos <j> ] 

It is then necessary that an extremizing trajectory satisfy 


3 L/ 


3v 


and similarly for y , <$ ,r, $ and 6 ; and also. 


3L/ 


3a 


= 0 


9L/ 

96 


= 0 


along with certain end conditions and transversality conditions which 
are not pertinent to this discussion. The latter two equations yield 


tan a = 


(tan 2 6 + 1 ) Ay 
v + *6 (tan 8 /cos Y) 


tan 6 = 


2 

(tan a + DA5 

( Ay tan a + v Av) cos y 


It would now be desirable to solve these equations for and 6 
so as to eliminate a and 8 from the equations of motion and the 
Euler-Lagrange equations. However, this is precisely the situation 
mentioned earlier; namely, the desired solutions cannot be readily 
obtained. 
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Now, the same original problem can be considered in a plumbline 

• i 

coordinate system with state variables x, y, z, x = 1, y = in, z = n, 
The relations among the various state variables in this and the 
previous system are given by: 


x = r cos (j) cos 0 

y = r cos <p sin 0 

z a r sin <p 

1 = v[sin y cos <p cos 0 - sin <f> cos 0 cos y cos 5 + 

sin 0 cos Y sin 6 ] 

m = v[sin Y cos <P sin 0 - cos 0 cos Y sin 5 - 
sin sin 0 cos Y cos 6 ] 

n = v[sin Y sin $ + cos 4> cos y cos 6 ] 


and, inversely, 


( 2 2 2 x 1/2 
v = ( 1 + m + n ) 


Y = 


6 = 


. -1 
sm 


xl + ym + zn 


cos 


-1 


(x 2 ty 2 + z 2 ) 1/2 (l 2 + m 2 + n 2 ) 1/2 

(x 2 +y 2 +z 2 )n - z(xl+ym+zn ) 
(x 2 +y 2 ) 1/2 ((xm-yl) 2 +(xn-zl) 2 +(yn-zm) 2 ) 1/2 


= (x 2 ♦ y 2 + z 2 ) 1 / 2 


. -1 
sin 


( x 2 + 7 V 77 k 


0 = tan* 1 Cy/x] 
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The equations of motion in the present coordinate system are: 


-F 


U x 

T~ „2 “ 2 ,3 7? 


M sin X cos X (x + y + z ) 

p y 


m = 




M cos X cos X (x 2 + y 2 + z 2)3/2 

p y 


n = F _ u z 

M sin X (x 2 + y 2 + 


x = X 
y = m 
z = n. 


Then, for 


L = a. 


-F 


yx 


M sin X cos X (x 2 + y 2 + z 2 ) 3/2 

— P y 


+ a 


m 


J2L 


M cos Xp cos Xy (x 2 + y 2 + z 2 ) 3 ^ 2 


+ a 


yz 


M sin X y (x 2 + y 2 + z 2 ) 3/2 _ 


+ a x l + a y m + a z n 


the Euler-Lagrange equations are: 
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cn = - o 

■L v 


a - . a 


ya, 


3 u x(x o 1 +y a^+z a n ) 


, 2 2 2 .3/2 

(x +y +z ) 


, 2 2 2.5/2 

(x +y +z ) 


pa 


3 yy(x 0l +y om+z on) 


, 2 2 2.3/2 

(x +y +z ) 


, 22 .5/2 

(x +y +z ) 


ya 


3U z(x a ^+y ° m +z ^ n ) 


, 2 2 2.3/2 

(x +y tz ) 


(xVtz 2 ) 5/2 


X_ cos X + o sin X cos X„ = 0 

P y m p y 


X sin X - a cos X sin X + 
p y m p y 

a cos X = 0 

n y 

along with the equations of motion. The latter two equations above 
are analogous to those obtained for the other coordinate system from 
3 L/3a =0 and 3 L/3B = 0. They give: 

tanX p = -°l /0 m , tan X y = o n /(o x 2 + o m 2 ) 1/2 

Substituting into the equations of motion gives: 


Aii 

*3 X 


a cos 


3 L 


3 in 
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F a. 


Wx 


M( a l 2+ a m 2+ a n 2)1/2 (x 2 +y 2 +z 2 ) 3/2 


F o 


m 


J2L 


M( a^2 + ° in 2 + ° n 2 )^ 2 (x 2 +y 2 +z 2 ) 3 ^ 2 


F o_ 


w z 


M ( a, 2 + o ra 2 + o n 2 ) 1/2 (x 2 + y 2 + z 2 ) 3/2 


X = 1 


y = m 


z = n 


Thus, the latter coordinate system provides an example of a 
situation in which the desirable solution for, and elimination of, the 
control angles is readily obtained. If a transformation from the 
latter coordinate system to the former coordinate system can now be 
effected, a similar set of solutions for a and ® should be immed- 
iately obtainable. 

Suppose, then, that the immediate example is set aside momentar- 
ily in order that a transformation technique can be discussed. This 
technique can then, hopefully, be applied to the example. Suppose a 
function L is given in a certain coordinate system with variables x^^ 

(i = 1, 2, ... , n) and Lagrange multipliers A. (i = 1, 2, , t . , n) 

as: 


L = x x , a 2 % • • • 9 ' 

Suppose, further, that a function L is given in another coordinate 
system with variables h- (i = 1, 2, ... , n) and Lagrange multi- 
pliers a. (i = 1, 2, ... , n), L = ^g^(t, h 1# h 2 » » h n ). 
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L leads to the Euler equations: 

3 f i 

X j = ~ X*^ /3 x_. (3 “ 1$ ••• ) n) 

and L to: 

. 3gi 

a ^ = - ck / 3hj (j = 1, ... , n) 

(where in both cases, the summation convention applies to the 
repeated sub-script i). 

Now, suppose that a transformation relationship is given for the 
variables in the two systems, say: 

h i = h i (x l* x 2» » x n } 

or, inversely: 

X| = x^(h^, h2i • •• 1 ) 


Then, for 


x i 



and 





f^(t , x^, . . . 



the function L = ^.f^ transforms into 


a. 

D 



( x 1 , . . , 


» ^ 


•i (t » X J 


x n } 
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The corresponding inverse multiplier transformation is: 



It can be demonstrated that, under these transformations, the 
Euler equations transform accordingly. Thus, if the functions 
gi(t, hp h 2 , ... , h n ) are properly related to the functions 
f.(t, x^, x 2 , ... , x n ) by the above transformation, the Lagrange 
multipliers are transformed as shown. 

In the example under consideration, all functions and transfor- 
mation relationships are properly defined for this application. Thus, 

3 v 3 v 95 

o, = A, V / 31 + A ■ / 31 + A / 31 

IV y o 

+ A 3r / 31 + A, 3(f y 31 + A a 30 / 3 1 
r <p o 

with corresponding equations for o , o a a and a where v, 

m n x y z 

are now expressed in terms of x, y, z, 1, m, n by 
the transformation equations. 

NOW, rhe various partial derivatives involved in preceding 

equations must be found in order to give the transformation equations 
explicitly. These partials can be found by differentiating the trans- 
formation equations with respect to the appropriate state variables and 
then solving the resulting algebraic system in the partial derivatives. 
Following this rather lengthy process will finally yield: 
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V 

1 



6 

1 


sin 9 cos Y sin 5 + cos 4> cos 9 sin Y - sin 4> cos 9 cos Y cos 5 

l/v[cos 4>cos 9 cos Y + sin <J> cos 9 s i n Y cos 

sin 9 sinY sin 5 1 

- ■ ■ ^ [sin 9 cos 5 + sin ♦ cos 9 sin 5 1 
v cos Y 

cos 4> sin 9 sin Y - cos 9 C os Y sin 5 - sin <P sin 9 C os Y cos 5 


Y = l/v[cos 9 sin Ysin 5 + cos 4> sin 9 C os Y + 
m 

sin <t> sin 9 sin Y cos 5 ] 

6 = - . ... 2: , [sin <l>sin 9 sin 5 - cos 9 C os 5 ] 

m v cos Y 

v = sin <P sin Y + cos <J> cos Y cos 5 
n 

Y = l/v[sin <f> cos Y - cos sin Y cos 5 ] 

5 = „ [cos $ sin 5 ] 

n v cos Y 

Further computation gives: 

2 2 2 22 2 2.2 2 
a l + °m + = x v + > Y / y + * 6 /v COS Y 

• • 

Now* the Euler equations for the plumbline system for 1, m, and 
n (in the form with tan X and tan X eliminated) may be transformed 

p y 

into the spherical system to yield an algebraic system of three 

« • • 

equations in v, Y * and 5 « This system may be solved algebraically 

• i • 

for v, Y , and 5 to gives 
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V = 


F Xy 

M( X 2 + X 2 y /v 2 + X 2 6 /v 2 cos 2 y ) 1/2 


g sm Y 


Y = 


F Xy + 

Mv( X 2 + X 2 ^ /v 2 + X 2 5 /v 2 cos 2 y ) 1 ‘^ 2 


(f- t) 


cos Y 


F A5 


Mv cos Y ( * 2 + * 2 $/v 2 + /v 2 cos 2 Y )^ 2 

v cos Ysin ^ 


r cot 4> 


These represent the Euler equations for the problem as expressed in the 
spherical coordinate system and with the control angles and 3 elim- 
inated. From them, or from direct transformation, 

tan a = ^Y /v ^ v 

tan S s X6 / v A ^cos Y 

This illustrates the transformation of Lagrange multipliers from 
one coordinate system to another and the accompanying transformation 
of the Euler equations. A particular value of such a transformation 
is that of obtaining expressions for the control angles in terms of 
the Lagrange multipliers in a coordinate system in which such ex- 
pressions could not readily be obtained in a direct fashion, as noted 
in the above. Other motivations for carrying out such a transformation 
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may also exist# For example, initial Lagrange multiplier values for 
one system may be used to give those for another. The simplicity and 
workability of such a transformation aid in making it a very valuable 
tool in trajectory optimization problems of the type illustrated. 
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Stability Criteria for n-th Order, Homogeneous 

f 

Linear Differential Equations 

E. F. Infante^ 

Center For Dynamical Systems, Brown University 





1. Introduction 

This note is concerned with the homogeneous differential equation 

x (n) + p.CtJx^* 1 ^ + ... + p n (t)x + P (t)x = 0, (1.1) 

1 n-1 n 

where the p^(t) are real continuous functions. It is desired to determine 
appropriate criteria for the stability of the origin, criteria dependent on 
the behavior of the functions p^(t) but not of their derivatives. 

This problem has been previously studied by Starzinski [1,2,3] for 
particular forms of this equation up to the fourth order, and by Razumichin 
[4] for the general matrix equation x = A(t)x. The approach of these 
authors has been to use the direct method of Liapunov, using a constant quad- 
ratic Liapunov function V(x) = x'3x which is generated by determining the 
n(n+l)/2 constant elements of the symmetric matrix B. The determination 
of all these elements requires very heavy algebraic computations, computa- 
tions which are completely unreasonable for n > 2. Recently, Ghizzetti 
[5,6] has obtained simple stability criteria for (1.1) by using some appro- 
priate maj oration formulae for all the integrals of this equation. The 


This research was supported by the National Aeronautics and Space 
Administration under Grant No. NAS8-11264. 
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particularly attractive aspect of these criteria is that they depend on only 
n constant parameters which locate a family of hyperellipsoids in the n- 
dimensional space of the p^(t). If the curve parametrically represented 
by the p^(t) is entirely contained within one of these ^hyperellipsoids , 
then (1.1) is asymptotically stable. 

In §2 of this note the second method of Liapunov is used to ob- 
tain stability criteria for (1.1) that depend on only n parameters which 
determine a family of elliptic paraboloids in the n-dimensional space p^(t). 
It can be shown that these elliptic paraboloids completely contain the hyper- 
ellipsoids of Ghizzetti, In § 3 a practical technique for the application of 
the stability criteria obtained is discussed and is applied in the last sec- 
tion to two examples. The stability conditions presented in this note are 
not necessary. Indeed, they are probably not the best possible conditions 
obtainable from a quadratic Liapunov function. The technique presented in 
this note was devised with particular emphasis on ease of computability of 
some simple criteria. 


2. Stability Criteria 

Consider Eq. (1.1) rewritten in state-space coordinates as 


x 

X 


1 

2 


x 

X 


2 

3 


x 


n-1 


= x 


n 


x 

n 


-P n <t)Xi 


... -p.(t)x . 
l n 


( 2 . 1 ) 
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It is assumed that the p^(t), real continuous functions of time, satisfy 
the Routh-Hurwitz inequalities [7]. Let the n real numbers a^, assumed 
to satisfy the Routh-Hurwitz inequalities, be associated to (2.1), which is 
rewritten as 


X 1 r X 2 


( 2 . 2 ) 


X = -(p (t) - o )x "... -(p/t) - a )x - ax - ... - ax. 
n n nl 1 lnnl In 


For economy of notation, (2.2) is rewritten as 


where 


x = Ax - U(t)x, 


(2.3) 


A = 


0 

1 

0 

. 0 

0 

0 

0 

1 

.. o 

• • 

0 

0 

0 

0 

9 • 

. 0 

1 

-a 

-a 

”01 - . 


-a_ 

i n 

n-1 

n-2 

2 

1 


U(t) 


0 \ 


/ n 

\- 


n 



n n-l 



# 

i » 

u = 

• 

u' ! 


1 n l / 


I u \ 
U \ 




n+l-i 

i 

u 


(2.4) 


and where n- = p.(t) - a.. 

i i l 


9b 


For the determination of the stability of the origin of (2,3), 

consider the Liapunov function V(x) = x r Bx, B* = B = (3^), 6^ = constant. 

Let b denote the n-th column of the matrix B, and 
n 


b 

n 





(2.5) 


The derivative V of the Liapunov function v in terms of (2.3) is given 
by 


V = x 1 ( A' B + BA)x - x * ( U 1 ( t )B + BU(t))x, 


( 2 . 6 ) 


or 


-V = x 1 Cx + x'(ub f t b u T )x, (2.7) 

n n 

where A'B + BA =-C. If it were possible to determine a matrix B, positive 
definite, such that -V is positive definite for all t 2 0, then asymptotic 
stability of the origin of (2.1) will have been determined by the well known 
theorem of Liapunov [8]. For this purpose, consider the following simple 
lemma : 

Lemma 2.1: Given the constant matrix A , defined by (2.4), for 

any constant positive semidefinite diagonal matrix C ? 0 the equation 
A T B + BA = -C has a unique solution B, and B is positive definite . 


94 



Proof. The matrix B, obviously symmetric, is unique since all the eigen* 


values of A have negative real parts. Now, let V(x^) = x q Bx o < 0 for 
some x q 4 0, and define as the trajectory of x = Ax issuing from 

x q at t = 0. Along 6 we then have V(x) s V(x q ) < ^ ut & approaches 
the origin and V(0) = 0. Hence V(x) > 0. Similarly, let V(x^) = °» 
x^ i 0, and 6^ the trajectory emanating from x^ at t Q . Since this 
trajectory approaches the origin, it must lie on the manifold x T Cx = 0. 

But this is clearly impossible with C diagonal and A in the form (2.4). 
Hence B is positive definite. 

Hence, let the matrix B be generated by the diagonal matrix 


C = 



( 2 . 8 ) 


u it 

where C and C are constant nonsingular positive definite diagonal 
square matrices, and where the zero element in the diagonal is located in the 
i,i position. On the basis of the above lemma V(x) = x'Bx will be posi- 
tive definite. In this case, Eq. (2.7) then becomes 



lc u | 


. u, u 1 u u T 

lu D + b u 
n n 

U u e .+n . .b u 

ni n+l-i n 

u, l' , u* V , 
u b + b u 
n n 

-V = x' 

0 

X + x' 

, u’ u f 

n ■. .b +8 .u 
n+l-i n xu 

28 .n . 

m n+l-i 

0 * II 

n ^ .b +8 .U 
n+l-i n ni 


1 C l l 


t JLu' J t u* 

\u b t b u 
‘ n n 

z z 

u 8 .+n -b 

ni n+l-i n 

Z.Z' £' 1 

u b + b u 

n n / 


(2.9) 
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Assume 3 ni >0 (it is always possible to find a 3^ > 0, namely 3 nn ) 
and consider the regular transformation x = Sy, 


S 



0 

1 

0 



( 2 . 10 ) 


where the unit element is in the i,i position and the I are unit matrices 
of appropriate dimensions. If this transformation is applied to Eq. (2.9), 
one obtains 


-V = y' 


y + y T 


„ .U , u 

3 .u -n , .b 
ni n+l-i n 


0 u 1 , u 1 

3 ,u -n .b 
ni n+l-i n 


26 niVl-i 


l l 

O' g .U -n , -b 
ni n+l-i n 


„ V .V 

g .u -n .b 
ni n+l-i n 


0 

( 2 . 11 ) 


or 


u' v u' 

-V = y' g .u - n . .b 
J ni n+l-i n 


g .u U -n .. -b U 

ru n+l-i n 


2 ^ni\+l-i 

3 .i/-n . .b*" 

ni n+l-x n 


0 \ 

o 1 Z ' 

3 .u - n , .b 

ni n+l-x n 


0 


It now becomes necessary to determine under what conditions (2.12) 
is positive definite. For this purpose, consider the second transformation, 
y = Tz, 


I 



T = 


O' 1 


O’ 




O' 



(2.13) 


where the unit element is in the i,i position, the I are unit matrices 

of appropriate dimensions and v u = S .u U - n , .b U , v^sg.u^-n . .b^. 

ni n+l-i n’ ni n+l-i n 

This transformation is obviously regular and when applied to Eq. (2.12) yields 


where 



(2.14) 


-1 

w = 2g .n . - (S .u u - n . .b u )’c u (8 .u u - n .b u ) + 
p m n+l-i ni n+l-i n ni n+l-i n 


, or 1 o o (2.15) 

(8 .u* - n .b x, )’C x ’ (6 .u fc - n . .b ) 

ni n+l-i n ni n+l-i n 


Since (2.14) is diagonal, it can then be concluded that V will be negative 
definite if u> t 6 > 0. 

On the basis of what has been said above, it is then possible to 

state: 


Theorem 2.1 : Given the homogeneous differential equation 
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x'*‘ # + Pl (t)x'“ _1 ' + ... + p _(t)x + p (t)x = 0 , (2. 

1 n-1 n 


with p.(t) real continuous functions for t i 0, associate with this 
equation the n real constants cu , . . . ,a satisfying the Routh- Hurwitz 


inequalities ) and define n ^ = p^(t) - cu. Let the matrix B = (B^) he 
the solution of the matrix equation 


A f B + BA 


( 2.17 


u # 

where C , C are constant, positive definite diagonal matrices, and the 


zero element in the diagonal 


m the 


Let denote the n-th column of B and define 


n n+1-1 | 

b = 8 . , u = n * . ' 

n ni n+l-i 





Then , if for any <5 > 0 and any i = l,...,n 


26 niVl-i 


<6 .u u - n ,b u )'C u 1 ( B .u“ - n -b u ) -t 
ni n+l-i n m n+l-i n 


( 2 . 20 ) 


(S .u 

ni 


,L,„i £ 

n+l-i n ni 


n .b*) > 6 
n+l-i n 


for all t > 0 , the null solution of (2,16) is asymptotically stable . 

This theorem is not as general as it would have been possible 

to state, yet it is still too general for practical applications because 

u j£ 

of the generality of the matrices C and C * Before restricting the 
theorem, it is desirable to make some remarks concerning the results so 
far obtained. 

First of all we wish to point out that Eq. (2.20) represents, 
in the parameter space of the n's, an elliptic paraboloid. This can be 
easily seen by introducing the transformation of coordinates for the para- 
meter space given by 



f W1 ' 




/ 6 .1 
ni 

-b U 

n 

0 \ 


i u \ 

V = 

Y n + l-i 

= 

Y n+l-i 


0* 

1 

0 ! 


n n+l-i 


l il 


l ^ / 


1 O’ 
* 

-b* 

n 

#1 ii 


l / 


This transformation is obviously regular if 6 n ^ > which as was pre- 

viously pointed out, is no restriction. In the new coordinates, Eq. (2.20) 
becomes 
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-1 


i cr 


23 -Y , . - v' 
ni 'n+l-i 


! v > e . 


( 2 . 22 ) 


.‘- 1 I 


This is evidently the equation of an elliptic paraboloid. If 8 ^ > 0 , as 
assumed, the domain defined in the parameter space by (2.22), hance by (2.20). 
is nonempty. 

U 

Secondly, it is evident that, for any C and C satisfying 
the conditions of Theorem 2 . 1 , the domain of the r) parameter space defined 
by any of the (2.20) is strictly contained within the domain where the p^(t) 
satisfy the Routh-Hurwitz inequalities. On the other hand, it is easily shown 
that every point of the domain of the parameter space where the p^(t) satis- 
fy the Routh-Hurwitz inequalities is contained in at least, one of the domains 
defined by ( 2 . 20 ). To prove this, let p^(t) n p i = c ° nst ants. Since the 

satisfy the Routh-Hurwitz inequalities, it is possible to select the n 
numbers ou, themselves satisfying these inequalities, and such that 
n ^ = P - a , . □ e > 0 for some j and p , . -a - . = 0 for 

all i i j. Under these conditions Eq. (2.20) reduces to 


23 .n .-n J . 1 .b u ’c u b u n . . - n , .b £ 'c £ b £ n . > 6. (2.23) 
n3 n+l-3 n+l-3 n n n+1-3 n+l-3 n n n+1-3 


But for any e > 0 sufficiently small, a 6 > 0 can be found such that ( 2 . 23 ) 
is satisfied. Hence the remark. 

Finally, it is noted that the continuity condition imposed by 
Theorem 2.1 on the p^(t) imply that Eq. ( 2 . 16 ) does not have a finite 
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escape tine. It is therefore possible on the basis of this remark and 
the two previous ones to state : 

Corollary 2.1: Given the differential equation (2.16) with 

p.(t) real continuous functions for t > 0, if there exist a t > 0 such 
that for all t > t (2.20) is satisfied for some 6 > 0 and some i = l,...,n, 
then the null solution of (2.16) is_ asymptotically stable . 

Corollary 2.2: If, in Eq . (2.16), the pAt) are real continuous 

functions for t > 0 and lim p^(t) = » where the p^ satisfy the 

t -*■ °° 

Routh-Hurwitz inequalities, then the null solution of (2.16) asymptotically 
stable . 

This last corollary is very well known [7], and can be traced 
directly to Liapunov. 

3. Application of Stability Criteria 

u 

The positive definite diagonal matrices C and C have not 
been so far specified. The first step in the application of the stability 
criteria obtained to a specific example is the selection of these two mat- 
rices, from which the matrix B is obtained as the solution of the equation 
A T B + BA = -C. Algorithms for the solution of this matrix equation are 
available. A particularly simple one has been recently given by Smith [9] 
in the case matrix A has the form (2.18). 

It is particularly convenient, to obtain algebraically simple forms 
for 3, to select the matrices C and C to oe composed ui liucar com- 
binations of matrices of the form 

C^ = 2 diag (p, 0 , . . . ,0) (3.1) 
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(3.2) 


C k = 2 diag (0,..., g-,0,...,0), k i 1, 

n 

where y is the Hurwitz determinant [7] of the a: 



a l 

a 3 

a 5 

• 



a 2n-l 


l 

a 2 

a 4 

• 



a 2n-2 


0 

a l 

a 3 

4 



a 2n-3 

y = 

0 

1 

a 2 

• 



a 2n-4 


• 

• 

• 




• 


• 

• 

• 


• 


• 


0 

0 

0 

• 

• 

• 

a 

n 


The matrix equation A'B^ t B^A = -C^, where A is given by (2.18) 

can be rapidly solved for when is of the suggested form. The 

matrices obtained in this manner for n = 2,3 are shown in Table 1. 

Ingwerson [10] previously published these matrices for n = 2,3,4. If 
u 

C and C are obtained, as suggested, by linear combinations of the C^, 
then the matrix B will be the corresponding linear combination of the B . 


Table I 


n = 2 


B, = 


a l + a 2 “l ! 


' 2a l a 2 0 


c = 


“1 

5 1 

[ 0 0 , 



P 

O 


0 0 1 

i 

CN 

1 



II 


C 2 = 


1 

| o 1, 


0 2a 


n = 3 



1 a 2 (a i a 2 _a 3 )+a i a 3 

2 

a i°2 

“iV 0 ^ 

i i 

1 2a 3 (a l° t 2 _a 3 ) 

0 

0 

3 1 = 

2 

a i a 2 

3 

“l + “3 

2 

a i 

i< 

H 

O 

r> 

0 

0 

0 


la, a„-a„ 

2 

o. 

i 

1 1 

l 0 

0 

0/ 


B 2 = 


1 a, a_ 


0 1 



1 ° 

0 

0 

1 3 

3 

l 



f 



“3 

2 

a^+a 

a i 

* 

C 2 = 

0 

2 ^ a i°‘2” a 3^ 

0 



i 



i 0 

0 

0 

n 

« 

a. 

i 1 





/ 



i°l 

a 2 a 3 




0 

0 1 

B 3 = 

a 2 a 3 

a l a 3 +a 2 

“3 

C 3 = 

0 

0 

0 


1 ° 

“3 

a 2 j 


U 

0 

2(a l a 2“ a 3 ) / 


103 



4. Two Examples 

In this section, the stability criteria obtained is applied to 
two simple but illustrative example problems. 

As a first example, consider the second order equation 

x + px + q(t)x = 0 (4.1) 


or 


x 

X 


1 

2 


x 


2 


-q(t)x 1 


p*. 


(4.2) 


where, p > 0 is a constant and o <q^+ £ < q(t) * q^ - £, for £ > 0. It is 

desired to determine conditions on q^, q^ and p that guarantee the 

asymptotic stability of the null solution of (4.2). This same problem has 

been treated by Ghizzetti [5], with whom we wish to compare our results. 

In the case of a second order equation, inspection of the matrices 

B and B of table one indicates that, for 6 . > 0 one must select i = 2 
12 ru 

With this choice one immediately obtains 


2a l a 2 


c = 


o o ! 



l 


2a l a 2 


y 


(4.3) 


' 0 2l| 

1 

1 

1 a l) 

1 I 

! n 2| 


1 q - a 2 


j 


II 

3 


= 


i 6 22; 


i 1 ! 

1 

l n i; 


1 p - a / 
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upon which the stability equation given by (2.20) becomes 


2(p - a 1 ) - [(q - a 2 ) - a^p - ap] 2 ^“ “ ° 2 * “ a / p " a l )] - 6 > 0 

(4.4) 


or, letting 




and z(t) = , 

P 


4V 1 V 2^ 1 “ v i^ 


[z(t) 


V l ( l " 


v x )] d e > 


(4.5) 


To determine the appropriate values of and for this 

expression, let 


z 


1 



\> 2 + v (1 - v x > - 2 /v^d - v 1 ) 


z 


2 


q 2 

— = v 2 + V l (1 “ v p + 2 *' v 1 v 2 (1 ~ V P 

P 


(4.6) 


and to maximize the difference between z^ and z^ let = 1/2. Then 


z 


1 





z 


2 




v 2 + 



(4.7) 


Solving now tor v rrom the firsi. uf 


v 


2 




(4.8) 


is obtained. With these two particular values of and v 2 (4.7) yields 
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z i + 2 ' / T* 


Z 1 + 




Hence, if 0 < z^ t £ < z(t) < £ for some £ > 0, an c > 0 
can be found such that Eq. (4.5) is satisfied. Therefore, Eq. (4.1) is 
asymptotically stable if, for some £ > 0, 


0 < q-L + 5 < q(t) < q 2 - c 


(4.10) 


and 


q 2 

2 

P 




q l 

2 

P 



(4.11) 


This result is represented in graphical form in Figure 1: if 

— 2 - is strictly internal to the domain A of the parameter space 
P 

2 2 

q^/p vs. q 2 /p , then Eq. (4.1) is asymptotically stable. The domain 
A obtained by Ghizzetti [5] is shown also. 

As a second example, consider the differential equation 


‘X t pX + x + r(t)x □ 0 , (4.11) 

where p > 0 is a constant and 0<£-r(t)-r 2 ~£ for some £ > 0. It 

is desired to determine conditions on r^ to guarantee the asymptotic sta- 
bility of the null solution of this equation. This equation has been 
studied by Starzinski [3], who generated a constant Liapunov function by 
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determining, through a very laborious process, appropriate values for all 

six elements of the 3x3 3 matrix. 

Inspection of the third order matrices of Table 1 indicates 

that, for 3 . > 0 one must select either i = 2 or i = 3. Let i = 3 
ni 

upon which the stability equation (2.20) becomes 

2633VP - a - [3 3 3(r(t) - - S 31 (p - °0] C 

(4.12) 

- C6 33 (l - a 2 ) - e 32 (p - a 1 )] 2 C £ * 6 > 0. 


Since i = 3, let C = where and are the two matrices 

shown in Table 1, and A > 0. From Table 1, then 


S 31 = °i 0l 2 " a 3» B 32 a 1 + Xa l’ B 33 


A + ou 


(4.13) 


-1 




-1 


2a 3^ a l a 2’ a 3^ 




2A(a^a2~a 3 ) 


are immediately obtained. Equation (4.12) can be therefore rewritten as 


4ot 3 (A + a 1 )(a 1 a 2 ”a 3 )(p - a^) - [(A + a 1 )(r(t) - c^) - ( a ]_ a 2" a 3^P “ + 

a (4 * 14 > 

- — C(X + a 1 >(l - a 2 > - (a 2 + Xap(p - ap] 2 > e > 0. 


The second quadratic term vanishes if 


(1 - op = opp - ap. 


(4.15) 
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Furthermore, (4.14) can be satisfied as r(t) becomes very small only if 


A + a. 


P - a = a 


1 3 a i a 2 " a 3 


(4.16) 


Assuming these two conditions, Eq. (4.14) yields 


0 < 5 s r(t) i 4 c* 3 - £ , 


(4.17) 


where £ 0 as e -*• 0. Equations (4.15) and (4.16) yield 


a 


3 



A + p 


P + 



4 + 4a 2 


(4.18) 


therefore , let 


a 2 = 1 - 

A) 2 

if 

0 < p < /2 


1 

a 2 " 2 

z 

if 

/2 < p 

(4.19) 


upon which one obtains that Eq. (4.11) is asymptotically stable if 


0 < C * r(t) S — — (p 2 + E_) - C if 0 < p £ /2 

(4.20) 

0 < C s r(t) < e if p > /2 

A + p ^ 


for some £ > 0 and X > q , 


since the a's obtained from Eq. (4.18) and 



(4.19) satisfy the Routh-Hurwitz inequalities. 

This same result would have been obtained if the stability 
Eq. (2.20) for i = 2 had been used. The stability conditions (4.20) are 
identical to those obtained by Starzinski [3]. 
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1 . Introduction . 

The purpose of this paper is to give a unified presenta- 
tion of Liapunov 1 s theory of stability that includes the classical 
Liapunov theorems on stability and instability as well as their 
more recent extensions. The idea being exploited here had its 
beginnings some time ago. It was, however, +he use made of this 
idea by Yoshizawa in [1] in his study of non autonomous differential 
equations and by Hale in [2] in his study of autonomous functional 
differential equations that caused the author to return to this 
subject and to adopt the general approach and point of view of this 
paper. This produces some new results for dynamical systems defined 
by ordinary differential equations which demonstrate the essential 
nature of a Liapunov function and which may be useful in applications. 
Of greater importance, however, is the possibility, as already in- 
dicated by Hale* s results for functional differential equations. 


*This research was support: e a in pari uy ihc National Aeronautics and 
Space Administration under Grant No. NGR^±^ under Contract 

No. NAS8-1126L, in part by the United States Air Force through the 
Air Force Office of Scientific Research under Grant No. AF-AFOSR-693-65 
and in part by the United States Army Research Office, Durham, under 
Contract No. DA-31-12L-ARO-D-270. 
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that these ideas can be extended to more general classes of dynam- 
ical systems. It is hoped, for instance, that it may be possible 
to do this for some special types of dynamical systems defined by 
partial differential equations. 

In section 2 we present some basic results for ordinary 
differential equations. Theorem 1 is a fundamental stability 
theorem for nonautonomous systems and is a modified version of 
Yoshizawa’ s Theorem 6 in [1], A simple example shows that the 
conclusion of this theorem is the best possible. However, when- 
ever the limit sets of solutions are known to have an invariance 
property then sharper results can be obtained. This "invariance 
principle" explains the title of this paper. It had Its origin for 
autonomous and periodic systems in [3] - [5], although we present 
here Improved versions of those results. Miller in [6] has estab- 
lished an invariance property for almost periodic systems and ob- 
tains thereby a similar stability theorem for almost periodic 
systems. Since little attention has been paid to theorems which 
make possible estimates of regions of attraction (regions of asymp- 
totic stability) for nonautonomous systems results of this type are 
included. Section 3 is devoted to a brief discussion of some of 
Hale’s recent results [2] for autonomous functional differential 
equations . 

2. Ordinary differential equations. 

Consider the system 
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X = f(t,x) 


( 1 ) 


n-i"! 

where x is an n-vector, f is a continuous function on R 
to R n and satisfies any one of the conditions guaranteeing unique- 
ness of solutions. For each x in R n we define |x| = 

2 2 — n 

(x, + . .. + x ) 2 , and for E a closed set in R we define 
v 1 n' > 

d(x,E) = Min {|x-y|: y in E}. Since we do not wish to confine our- 
selves to hounded solutions, we introduce the point at » and 
define d(x,») = |xj~^ . Thus when we write E* = E {j{°° }, we shall 
mean d(x,E*) = Min{d(x,E), d(x,«)}. If x(t) is a solution of 
(1), we say that x(t) approaches E as t -> » if d(x(t),E) -> 0 
as t -» ». If we can find such a set E, we have obtained in- 
formation about the asymptotic behavior of x(t) as t -> ». The 
best that we could hope to do is to find the smallest closed set 
£2 that x(t) approaches as t -»», This set £2 is called the 
positive limit set of x(t) and the points p in £2 are called 
the positive limit points of x(t). In exactly the same way one 
defines x(t) -> E as t -» -» , negative limit sets, and negative 
limit points. This is exactly G. D. Birkhof f f s concept of limit 
sets. A point p is a positive limit point of x(t) if and only 
if there is a sequence of times t^ approaching ® as n -* » and 
such that x (t n ) P as n -> « . In the above it may be that the 
maximal interval of definition of x(t) is [0 ,t) . This causes 
no difficulty since in the results to be presented here we need 
only with respect to time t replace «by t. We usually ignore 
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this possibility and speak as though our solutions are defined on 
[0,oo) or (- 00 , 00 ) . 

Let V(t,x) he a C 1 function on [0,°°) x R n to R, and 
let G he any set in R n . We shall say that V is a Liapunov 
function on G for equation (l) if V(t,x) ^ 0 and V(t,x) g 

-W(x) ^ 0 for all t > 0 and all x in G where W is 
continuous on R n to R and 

v = sv + s av f . . 

St i=l Sx. 1 

1 

We define (G is the closure of G) 

E = {x; W(x) = 0, x in G] . 

The following result is then a modified hut closely re- 
lated version of Yoshizawa* s Theorem 6 in [1], 

THEOREM 1. If V is a Liapunov function on G for equation (l), 
then each solution x(t) of (l) that remains in G for all 
t > t ^ 0 approaches E* = E U {«>} as t -» <» , provided one of 
the following conditions is satisfied: 

(i) For each p in G there is a neighborhood 1\| of 
p such that |f(tjx)| is hounded for all t > 0 and 
all x in N . 

(ii) W is and W is hounded from above or below 

along each solution which remains in G for all 

t > t 1 0 . 
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If E is bounded, then each solution of (1) that remains in G 
for t > t § 0 either approaches E or » as t -» <» . 

Tims this theorem explains precisely the nature of the 
information given by a Liapunov function. A Liapunov function 
relative ter a set G defines a set E which under the conditions 
of the theorem contains (locates) all the positive limit sets of 
solutions which for positive time remain in G. The problem in 
applying the result is to find "good" Liapunov functions. For 
instance, the zero function V = 0 is a Liapunov function for the 
whole space R n and condition (ii) is satisfied but gives no in- 
information since E = R n . It is trivial but useful for appli- 
cations to note that if and are Liapunov functions on G, 
then V = + V 2 is also a Liapunov function and E = E^ f! E^ . 

If E is smaller than either or , then V is a "better" 
Liapunov function than either V-]_ or V 2 and is always at least as 

"good" as either of the two. 

Condition (i) of Theorem 1 is essentially the one used 

by Yoshizawa. We now look at a simple example where condition (ii) 

is satisfied and condition (i) is not. The example also shows that 

the conclusion of the theorem is the best possible. Consider 

2 2 

x + p (t ) x + x = 0 where p(t) ^ 8 > 0 . Define 2V = x + y y 

2 2 

where y = x . Then V = -p(t)y ^ - 8y and V is a Liapunov 
function on R 2 . Now W = 8y 2 and W = 28yy = -26 (xy + p(t);y 2 ) £ 
-28xy. Since all solutions are evidently bounded for all t > 0, 
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condition (ii) is satisfied. Here E is the x-axis (y = 0) 
and for each solution x(t), y(t) = x(t) -* 0 as t -» « . Noting 
that the equation x + (2 + e"^)x + x = 0 has a solution 
x (t ) = 1 + e~"k , we see that this is the best possible result with- 
out further restrictions on p . 

In order to use Theorem 1 there must be some means of 
determining which solutions remain in G . The following corollary 9 
which is an obvious consequence of Theorem 1, gives one way of 

doing this and also provides for nonautonomous systems a method for 
estimating regions of attraction. 

Corollary 1. Assume that there exist continuous functions u(x) 
and v (x) on R n to R such that u(x) § V(t,x) g v(x) for all 
t S 0 . Define Q* = {x ; u(x) < tj] and let G + he a component 
of Q + . Let G denote the component of Q = (x ; v(x) < rj} 
containing G + . If V is a Liapunov function on G for (1) and 
the conditions of Theorem 1 are satisfied, then each solution of 
(1) starting in G + at any time t Q ^ 0 remains in G for all 
t > t Q and approaches E* as t . If G is bounded and 
E° = E n G C G + , then E° is an attractor and G + is in its 
region of attraction. 

In general we know that if x(t) is a solution of 
(l)--in fact, if x(t) is any continuous function on R to R n -- 
then its positive limit set is closed and connected. If x(t) is 
bounded, then its positive limit set is compact. There are, how- 
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ever, special classes of differential equations where the limit 
sets of solutions have an additional invariance property which 
makes possible a refinement of Theorem 1. The % first of these are 
the autonomous systems 

X = f(x) (3) 

The limit sets of solutions of (3) are invariant sets. If x(t) 
is defined on [0,°°) and if p is a positive limit point of x(t), 
then the points on the solution through p on its maximal inter- 
val of definition are positive limit points of x(t). If x(t) is 
hounded for t > 0 , then it is defined on [0,»), its positive 
limit set ft is compact, nonempty and solutions through points 
p of ft are defined on (-<»,») (i.e., ft is invariant). If 
the max imal domain of definition of x(t) for t > 0 is finite, 
then x(t) has no finite positive limit points; that is, if the 
maximal interval of definition of x(t) for t > 0 is [0,P), 
then x(t) -» oo as t -» g . As we have said before, we will always 
speak as though our solutions are defined on (-°°,°°) and it should 
be remembered that finite escape time is always a possibility unless 
there is, as for example in Corollary 2 below, some condition that 
rules it out. In Corollary 3 below, the solutions might well go to 
infinity in finite time. 

The invariance property of the limit sets of solutions 
of autonomous systems (3) now enables us to refine Theorem 1. 

Let V be a C 1 function on R n to R . If G is any arbitrary 
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set in R n , we say that V is a Liapunov function on G for 
equation (3) if V = (grad V)’ f does not change sign on G . 
Define E = { x ; V(x) = 0 , x in Gf j , where G is the 
closure of G . Let M be the largest invariant set in E . M 
will be a closed set. The fundamental stability theorem for 
autonomous systems is then the following: 

THEOREM 2. If V is a Liapunov function on G for (3), then 
each solution x(t) of (3) that remains in G for all t > 0 
(t < 0) approaches M* = M U {»} as t -» « (t -> -») . If M is 
bounded, then either x(t) -» M or x(t) -» co as t -» « (t -> -°°) . 

This one theorem contains all of the usual Liapunov like 
theorems on stability and instability of autonomous systems. Here 
however, there are no conditions of definiteness for V or V , 
and it is often possible to obtain stability information about a 
system with these more general types of Liapunov functions. The 
first corollary below is a stability result which for applications 
has been quite useful and the second illustrates how one obtains 
information on instability. Cetaev* s instability theorem is 
similarly an immediate consequence of Theorem 2 (see section 3)* 

COROLLARY 2. Let G be a component of Q = { x ; V(x) < t) } . 
Assume that G is bounded, V ^ 0 on G , and M° = M h G C G . 
Then M° is an attractor and G is in its region of attraction. 

If, in addition, V is constant on the boundary of M° , then 
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M° is a stable attractor. 

Note that if M° consists of a single point p , 
then p is asymptotically stable and G provides an estimate of 
its region of asymptotic stability. 

COROLLARY 3. Assume that relative to ( 3 ) that V V > 0 on G 
and on the boundary of G that V = 0 . Then each solution of 

(3) starting in G approaches » as t -> «> (or possibly in 

finite time). 

There are also some special classes of nonautonomous 
systems where the limit sets of solutions have an invariance 
property. The simplest of these are periodic systems (see [3]). 

x = f(t,x) , f(t + T,x) = f(t) for all t and x . (k) 

Here in order to avoid introducing the concept of a periodic 

approach of a solution of (4) to a set and the concept of a 

periodic limit point let us confine ourselves to solutions x(t) 
of (1) which are bounded for t > 0 . Let ft be the positive 
limit set of such a solution x(t), and let p be a point in ft . 
Then there is a solution of (4) starting at p which remains in 
ft for all t in (- 00 , 00 ) ; that is, if one starts at p at the 

proper time the solution remains in ft for all time. This is the 

sense now in which ft is an invariant set. Let V(t,x) be C 1 
on R x k 31 and periodic in * of period T . For an arbitrary 
set G of R n we say that V is a Liapunov function on G for 
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for the periodic system (4) if V does not change sign for all 
t and all x in G . Define E = { (t,x) ; V(t,x) = 0, x in G 
and let M be the union of all solutions x(t) of (4) with the 
property that (t,x(t)) is in E for all t . M could be called 
"the largest invariant set relative to E" . One then obtains the 
following version of Theorem 2 for periodic systems: 

THEOREM 3. If V is a Liapunov function on G for the periodic 
system (4), then each solution of (4) that is bounded and remains 
in G for all t > 0 (t < 0) approaches M as t -> oo (fc -> -») . 

In [6] Miller showed that the limit sets of solutions 
of almost periodic systems have a similar invariance property and 
from this he obtains a result quite like Theorem 3 for almost 
periodic systems. This then yields for periodic and almost periodic 
systems a whole chain of theorems on stability and instability 
quite similar to that for autonomous systems. For example, one has 

COROLLARY 4. Let Q+ = ( x; V(t,x) < rj, all t in [0,T] } , and 
let G + be a component of Q* . Let G be the component of 

= { x; V(t,x) < rj for some t in [0,T] } containing G + . If G 
is bounded, V ^ 0 for all t and all x in G , and if M° = 

M fl G C G + , then M° is an attractor and G + is in its region of 
attraction. If V(t,x) = cp(t) for all t and all x on the 
boundary of M° , then M° is a stable attractor. 

Our last example of an invariance principle for ordinary 
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differential equations is that due to Yoshizawa in [1] for "asymp- 
totically autonomous" systems. It is a consequence of Theorem 1 
and results by Markus and Opial (see [1] for references) on the 
limit sets of such systems. A system of the form 

X = F(x) + g(t,x) + h(t,x) (5) 

is said to be asymptotically autonomous if (i) g(t,x) -» 0 as 
t -> oo uniformly for x in an arbitrary compact set of R 11 , 

oo 

(ii) j |h(t,cp(t))| dt < » for all cp bounded and continuous 
o 

on [0,») to R . The combined results of Markus and Opial then 
state that the positive limit sets of solutions of (5) are in- 
variant sets of x - F(x) . Using this, Yoshizawa then Improved 
Theorem 1 for asymptotically autonomous systems. 

. It turns out to be useful, as we shall illustrate in a 
moment on the simplest possible example, in studying systems (1) 
which are not necessarily asymptotically autonomous to state the 
theorem in the following manner: 

THEOREM U. If, in addition to the conditions of Theorem 1, it is 
known that a solution x(t) of (l) remains in G for t > 0 
and is also a solution of an asymptotically autonomous system (5), 
then x(t) approaches M* = M U (oo) as t , where M is the 

largest invariant set of x = F(x) in E . 

It can happen that tne system (l) is itself asymptotically 
autonomous in which case the above theorem can be applied. However, 
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as the following example illustrates , the original system may not 

itself be asymptotically autonomous but it still may be possible 

to construct for each solution of (l) an asymptotically autonomous 

system (5) which it also satisfies. 

Consider again the example 

x = y (6) 

y = -x - p(t)y , 0 < 6 § p(t) S m 

for all t > 0 

Now we have the additional assumption that p(t) is bounded from 
above. Let (x(t), y(t)) be any solution of (6). As was argued 
previously below Theorem 1, all solutions are bounded and y(t) -» 0 
as t -» oo . Now (x(t), y(t)) satisfies x = y , y = 

-x - p(t)y(t), and this system is asymptotically autonomous to 
(*) x =i y , y = -x . With the same Liapunov function as before, 

E is the x-axis and the largest invariant set of (*) in E is the 
origin. Thus for (6) the origin is asymptotically stable in the 
large . 

3. Autonomous functional differential equation. 

Difference differential equations of the form 

x(t) = f (t,x(t),x(t-r)) , r > 0 (7) 

have been studied almost as long as ordinary differential equations 
and these as well as other types of systems are of the general form 


where x is in R n and x^ is the function defined on [-r,0] 
by x (t) = x(t+r), -r S t § 0. Thus x is the function that 
describes the past history of the system on the interval [t-r,t] 
and in order to consider it as an element in the space C of 
continuous functions all defined on the same interval [-r,0], x^_ 

is taken to he the function whose graph is the translation of the 
graph of x on the interval [t-r,t] to the interval [~r,0] . 

Since such equations have had a long history it seems surprising 
that it is only within the last 10 years or so that the geometric 
theory of ordinary differential equations has been successfully 
carried over to functional differential equations. Krasovskii [8] 
has demonstrated the effectiveness of a geometric approach in ex- 
tending the classical Liapunov theory, including the converse 
theorems, to functional differential equations. An account of other 
aspects of their theory which have yielded to this geometric approach 
can be found in the paper [9] hy Hale. What we wish to do here is 
to present Hale 1 s extension in [2] of the results of Section 2 of 
this paper to autonomous functional differential equations 

X = f(x t ) . (9) 

It is this extension that has had so far the greatest success in 
studying stability properties of the solutions of systems (9), and 
it is possible that this may lead to a similar theory for special 
classes of systems defined by partial differential equations. 

With r 1 0 the space C is the space of continuous 
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functions cp on [-r,0] to R n with |jcp|| = 

max { | cp (nr) | ; -r ^ t ^ 0} . Convergence in C is uniform conver- 
gence on [-r,0]. A function x defined on [-r,<») to R n is 
said to be a solution of ( 9 ) satisfying the initial condition cp 

at time t = 0 if there is an a > 0 such that x(t) = f(x ) 

for all t in [0,a) and x q = cp . Remember x q = cp means 
x(t) = cp(r), -r ^ t ^ 0. At t = 0, x is the right hand deriv- 
ative. The existence uniqueness theorems are quite similar to 
those for ordinary differential equations. If f is locally 
Lipschitzian on C, then for each cp in C there is one and only 
one solution of ( 9 ) and the solution depends continuously on cp . 

The solution can also be extended in C for t > 0 as long as it 
remains bounded. As in Section 2 ; we will always speak as though 
solutions are defined on [-r,oo). The space C is now the state 
space of ( 9 ) and through each point cp of C there is the motion 

or flow x starting at cp defined by the solution x(t) of (9) 

satisfying at time t = 0 the initial condition cp; x^_, 0 i t O, 
is a curve in C which starts at time t = 0 at cp. In analogy 
to Section 2 with C replacing R n , x^ replacing x(t), and 
||xj| replacing |x(t)|, we define the distance d(x^,E) of x^_ 

from a closed set E of C to be d(x ,E) = min {||x € E} . 
The positive limit set of x^_ is then defined in a manner completely 
analogous to Section 2. Because there are some Important differences 
we shall be satisfied here with restricting ourselves to motions 
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x t bounded for t > 0, One of the differences here is that in 
C closed and hounded sets are not always compact. Another is that 
although we have uniqueness of solutions in the future two motions 
starting from different initial conditions can come together in 
finite time t >0: after this they coincide for tit . flhe 

motions define semi-groups and not necessarily groups.) 

Hale in [2] has, however, shown that the positive limit 
sets ft of hounded motions are nonempty, compact, connected, 

invariant sets in C . Invariance here is in the sense that, if 
x^ is a motion starting at a point of ft, then there is an exten- 
sion onto (-oo, -r] such that x(t) is a solution of (9) for all 
t in (-oo^oo) and remains in ft for all t . With this 

result he is then able to obtain a. result which is similar to 
Corollary 1 of Section 2 # 

For cp e C let x (cp) denote the motion defined by (9) 
starting at cp . For V a continuous function on C to R define 

V and by 

V(cp) = lim i [V(x (<p))-V(cp)]. (10) 

T -> Of 

and 

Qg = (cp ; V (cp) < &} . 

THEOREM 5. If V is a Liapunov function* on G for (9) and x is 

u 

a trajectory of (9) which remains in G and is bounded for t > 0, 
then x , — > M as t . 

* As before, V is a Liapunov function on G, if V does not 
change sign on G. 
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Hale has also given the following more useful version 
of this result. 

COROLLARY 5. Define Q = £cp ; V(cp) < rj} and let G he or 

a component of . Assume that V is a Liapunov function on 
G for (9) and that either (i) G is bounded or (iii) |cp(0)| is 
bounded for cp in G . Then each trajectory starting in G 
approaches M as t -» 00 . 

The following is an extension of Cetaev T s instability 
theorem. This is a somewhat simplified version of Hale 1 s Theorem 4 
in [2], which should have stated ,f V(cp) > 0 on U when cp / 0 

and V( 0 ) = 0 " and at the end "... intersect the boundary of 

This is clear from his proof and is necessary since he 
wanted to generalize the usual statment of Cetaev* s theorem to in- 
clude the possibility that the equilibrium point be inside U as 
well as on its boundary. 

COROLLARY 6 . Let p € C be an equilibrium point of ( 9 ) contained 

in the closure of an open set U and let N be a neighborhood of 

p . Assume that (i) V is a Liapunov function on G = U (1 N, 

(ii) M fl G is either the empty set or p, (iii) V(cp) < r\ on G 

when cp £ p, and (iv) V(p) = r\ and V(cp) = r\ on that part of 

the boundary of G inside N. Then p is unstable. In fact, if 

N is a bounded neighborhood of p properly contained in N then 
o 

each trajectory starting at a point of G q = G H N q other than p 
leaves N in finite time. 
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Proof o By the conditions of the corollary and Theorem 6 each 


trajectory starting inside G q at a point other than p must 

either leave G^, approach its boundary or approach p . 

Conditions (i) and (iv) imply that it cannot reach or approach that 

part of the boundary of G q inside N q nor can it approach p 

as t w . Now (ii) states that there are no points of M on 

that part of the boundary of N q inside G . Hence each such 

trajectory must leave N in finite time. Since p is either in 

o 

the interior or on the boundary of G, each neighborhood of p 
contains such trajectories, and p is therefore unstable. 

In [2] it was shown that the equilibrium point cp = 0 of 

&(t) = ax^(t) + bx^(t-r) 

was unstable if a > 0 and (b| < j a| . Using the same Liapunov 
function and Theorem 5 we can show a bit more. With 

v(<p) = - + t / q> 6 (e)ae , 

4a -r 

v(x, ) = - | / x b (e)ae 

4a t-r 

and 

V(q>) = -4(9 6 (0) + 2 5 <p 3 (0)cp 3 (-r) + q> 6 (-r)) 

d 

which is nonpositive when |b| < | a| (negative definite with re- 
spect to cp(0) and cp(-r)); that is, V is a Liapunov function 
on C and E - [cp; qp (0) = qp(-r) = 0} . Therefore M is simply 
the null function qp = 0 . If a > 0, the region G = (cp; V(cp) < 0} 
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is nonempty, and no trajectory starting in G can have cp = 0 as 
a positive limit point nor can it leave G . Hence by Theorem 5 
each trajectory starting in G must be unbounded. Since qp = 0 
is a boundary point of G, it is unstable. It is also easily seen 
[2] that if a < 0 and |b| < j a| , then cp = 0 is asymptotically 
stable in the large. 

In [2] Hale has also extended this theory for systems 
with infinite lag (r = »), and in that same paper gives a number 
of significant examples of the applications of this theory. 
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An analytical solution of the Euler-Lagrange equations for the Lagrange 
multipliers for optimum coast trajectories (minimum fuel consumption) 
is obtained. Previous solutions have a singularity at zero eccentri- 
city. The present solution does not have this singularity, but there 
is a numerical difficulty due to a removable singularity at unit 
eccentricity. An approximate solution, accurate near ’unit eccentricity, 
is given. This solution reduces to the exact parabolic solution for 
unit eccentricity. 
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1. INTRODUCTION 


Optimization of the flight trajectory of a rocket powered space vehicle 
with the indirect method of calculus of variations requires the simul- 
taneous integration of the equations of motion (constraint equations) 
and the Euler-Lagrange equations for the Lagrange multipliers. Because 
of the difficulty in obtaining an analytical solution to this problem 
during powered portions of the trajectory, the integration must be 
performed numerically. For coasting portions (zero thrust) of the tra- 
jectory an analytical solution for the motion of the vehicle can be 
produced if it is assumed that the vehicle moves in a vacuum under the 
action of the gravitational field of a single central body (spherical 
earth). The motion is, of course, governed by the classical Kepler 
solution of the two body problem. However, the Lagrange multipliers are 
not constant during coasting portions of the trajectory and it is neces- 
sary to solve the Euler-Lagrange equations for the multipliers to deter- 
mine the optimal direction of thrust at the end of each coasting arc. 

It is the purpose of this paper to present a new analytical solution for 
the multipliers. This solution used with the appropriate form of Kepler’s 
solution gives a particularly convenient analytical form of solution for 
optimal coasting arcs. When implemented in a numerical routine for 
trajectory computation, the analytical solution not only reduces compu- 
tational time, but also eliminates the errors due to numerical integra- 
tion for long coasting arcs. In studies made of the earth orbit rendez- 
vous problem 1 , it was found that long coasting arcs were often a 
necessity to avoid severe payload penalties. In this case, the analytical 
solution for coast is of especial advantage. 

Analytical solutions for coast were presented by W. E. Miner^ in 1963 
and by S. A. Jurovics3. In the June, 1965* issue of the AIAA Journal , 

M. W. Eckenwiler^ presented a solution very similar to Miner's. These 
solutions have the disadvantage of a singularity at zero eccentricity. 
Since circular and near circular orbits are of major interest, this 
singularity is removed in the present solution. A numerical difficulty 
appears here for near unit eccentricity. However, an exact solution for 
parabolic orbits and, in addition, an approximate solution for near 
parabolic orbits that eliminates the numerical difficulty are given. 


2. MINIMUM FUEL TRAJECTORY 

The equations of motion for a rocket in the gravity field of a spherical, 
homogeneous earth are: 

TT a V 



( 1 ) 
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2. (Continued) 


The state variables are the position vector R, the velocity V, and the 
mass m. The control variable is the thrust vector T. The gravitational 
constant of the central boyd and the magnitude of F are denoted by y 
and r, respectively. The mass flov rate of the rocket is assumed to be 
constant, and given by 

m = -M (2) 

For minimum fuel consumption, it is required that 

- m dt = MINIMUM (3) 

Appending the equations of motion, Eqs. (l) and the mass flow require- 
ment Eq. (2) as constraints with the Lagrange multiplier techniques, 
the following integral must be a minimum 

* tf • — 2- — 

(-m + X • (V + -^ F - -) + 7 • (R - V) + o(m + M)]dt (4) 

$ = r 3 m 

J t 

o 

The Components of the vectors X and Y, and a are the Lagrange multipliers. 
The Euler-Lagrange equations for these Lagrange multipliers are: 

m 

X + Y * 0 


y . -L X + 2 m. (a • R)F = 0 
r* 5 r^ 



m 


where T and X are the magnitudes of T and X, respectively. The third of 
these equations shows that X is the same direction as the thrust. When 
the thrust is zero, the last equation implies that a is a constant. The 
following differential equation for X comes from the first two of Eqs. (5)* 

X — — X + (X • R)R (6) 

r -5 

This equation and the constraint equations have to he satisfied along 
the optimum trajectory. 
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3. ANALYTICAL SOLUTION FOR LAGRANGE MULTIPLIERS DURING COAST PERIOD 


During coast the plane of motion is fixed, and from Eq. (6) the differ- 
ential equation for the component of A normal to the plane of motion 
"becomes 



( 7 ) 


Let x and y represent the components of the Kepler solution in a carte- 
sian coordinate system with the x axis directed toward perigee and the 
y axis in the plane of motion. By comparing Eq. (7) with the equation 
of motion during coast, we see that the solution of Eq. (7) is 


X N = K x x + K 2 y (8) 

where and K 2 are the integration constants. It is particularly _ 
convenient to assume the solution of Eq. (6) for the projection of A 
on the plane of motion to have the form 

Ap = FR + GR (9) 

The F and G are, in general, functions of time and uniquely define Ap 
when the position and the velocity vectors do not coincide. 

Note that the form of Eq. (9) is similar to that assumed for the position 
vector in terms of the classical f and g series^. Here, F and G are not 
these series, hut, as shall be found, finite expressions involving 
quantities from the Kepler solution. 

Substituting Eq. (9) into Eq. (6), we have after manipulation 


•• • • 

Ap = 2L- (2F + 3^ G)R - GR (10) 

r 3 r r 3 

Eq. (9) is differentiated twice with respect to time. This yields: 

Xp = (F + ^HrG - + (2F — ^-G + G)R (ll) 

r 4 r>3 jO r-3 

By domparing Eqs. (10) and (ll), we see that: 

G = -2F 

P = ^ (3F + 2G) 
r J 


134 


( 12 ) 

(13) 



3* (Continued) 

From Eq. (12), we have 


G = K 3 - 2F [lb 

By substituting Eq. (lU) into Eq. (13) > a differential equation for F 
alone is obtained 


P + ^ (F - 2 K 3 ) = 0 


This equation has the solution 

F = K^X + K 5 Y + 2K 3 (16) 

where, again, x and y are the components of the Kepler sblution and K^, 
are integration constants* The function G is obtained by inte- 
grating Eq. (l4). 

G - |(K 3 - 2F)dt = - j(3K 3 + 2Kl|X+ 2K 5 Y)dt 

In order to carry out this integration, the eccentric anomaly of the 
Keplerian motion is used as the integration variable for elliptical 
orbits. 

G = -3K 3 t - 2Ki * | —( cos u - e)(l - e cos u)du - 2 K 5 J ^(1 - e 2 )^ 
sin u(l - e cos u)du 
The integration gives 

G = - 3 Kot + ±K U + P) - 3et] + P(r + P) + Kg (IT) 


-'3 2 H E L d 

where t = time 

a = semi -major axis 
n = mean motion 
E = total energy 
L = angular momentum 
P = semi-latus rectum 
e = eccentricity 
Kg = integration constant 

The G function has exactly the same iorm ±u± hyperbolic orbit Q y*hf>r\ the 
integration is performed by using hyperbolic anomaly as the integration 
variable. Similarly, for parabolic orbits, the G function becomes 

G = -3K 3 t + 1^ £[r 2 - P(r + P)3 + K 5 £ P'( r + P) + Kg (l 8 ) 
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3 . (Continued) 

Finally, Y can “be obtained from 

• M » 

7 = -T = -GR + LD + K 3 R - (K-jX + K 2 Y)k ( 19 ) 

with 

D = -K 5 i + K4 J 

where i, J and k are unit vectors of the perigee oriented coordinates. 


It. APPROXIMATE SOLUTION OF COAST FOR ORBITS WITH ECCENTRICITY VERY CLOSE 
TO UNITY 


The elliptical solution as found in the previous section is 
A n = K-jX + K 2 Y 
Ap = FR + GR 
F = Kl^X + K 5 Y + 2 K 3 

G = - 3 K 3 t + [I(r + P) - 3et] + K^>(r + P) + K6 


( 8 ) 

( 9 ) 

(16) 

(IT) 


All terms are well behaved numerically as eccentricity approaches unity 
except the term with coefficient in the G- function. The energy E and 
the squared bracket term vanish separately at unit eccentricity but their 
ratio is finite. To eliminate this numerical difficulty, let eccentri- 
city e = 1 - e with e << 1, then the eccentric anomaly u is also much 
less than 1 and can be expanded in a series of sin u, which is: 

u = sin u + ~ • i sin^u + i . 3 . • — sin^u + .... (20) 

2 J 2 U 5 

This series is truncated after the third term and substituted into 
Kepler* s equation: 


t = — (u - e sin u) (21> 

ri 

This yields 

t = —[sin u + ^sin^u + ^sin^u - (l - e)sin u] 

= rKe sin u + ^sin^u + ^sin^u] 
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1+. (Continued) 

Since 

y 

sin u = ^<1 


1 


e 2 ) 2 " 


Y 


P 


1 

e 2 (2 



and 

i iii 

n = (\) 2 = JWl - e 2 ) 2 = £- c 2 (2 -e) 2 

p 2 


Kepler's equation becomes 



+ 


1 y£ 

2 P 


±_ Y U , 

+ £ lo p3 ] 


( 22 ) 


The second and higher order terms in e have been dropped in the last 
expression for time. These terms will also be neglected in the deriva- 
tion hereafter. 


Now 


Since 


with 


thus 


Y 2 = a^(l - e 2 ) sin 2 u = aP(l + cos u)(l - cosuu) 

2 

= (1 + cos u)a[l - (e + e)cos u] 

= (1 + cos u)(r - ae cos u) 

cos u = i (1 - £) = (1 + e) (1 - £) - 1 + e(l - 2£) 
e a a a 

« - |U ♦ §> 

( 2 r - P) ♦ 2 e[ (2r - P) - j£] 


Substituting Eq. (23) into Eq. (22), yields 

i - I- (x : T) * |[(r + P) + k£] 

3L 5 * 


(23) 


(2k) 
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(Continued) 


With this expression for time, the term with coefficient in the 
G- function is 

E {¥ r + p ) - (1 - £ ) I Ur + P) + |<r + P + !^)]j 
= |k u X { r 2 « p(r + p ) - i.[P(r + P) + Ifr 2 ]} 

Finally we have 

G = -3K 3 t + |K 4 I |r 2 - P(r + P) - J[P(r + P) + 1+r 2 ]} + (25) 

K 5 ^P(r + P) + Kg 

When e = 0, this becomes the parabolic solution given by Eq. (18). 


5. DISCUSSION 


The singularity at zero eccentricity that exists in previous analytical 
coast solutions does not appear in the solution presented here. The 
form of solution taken for the projection: of X onto the plane of motion 
eliminated this singularity and also gave a particularly simple form for 
the resulting solution. The absence of the singularity is of assistance 
for study of optimal trajectories for boost to circular and near cir- 
cular orbits. With the improved accuracy due to elimination of integra- 
tion and round off errors, the analytical solution presented is also 
useful for optimizing parking orbits. 

The numerical difficulty for very nearly parabolic orbits can be handled 
by employing the approximate solution presented in Section h. It is of 
interest to note, however, that in actual numerical work carried out, 
the exact parabolic solution gives satisfactory results in that region. 
The average of the significant figures of agreement of the Lagrange 
multipliers obtained from the analytical solution with the multipliers 
calculated by direct numerical integration is shown in the figure below. 
The comparison was made 100 seconds after the beginning of coast. Coast 
was initiated 1000 seconds after perigee. The values were calculated 
carrying eight digits. In the ranges of eccentricity away from unity, 
the difference in analytical and numerical values is due to error in 
the numerical integration. In these regions, there is no numerical 
difficulty with the analytical solution and it gives precise values. 

The results in the neighborhood of eccentricity e * 1 *10“^ can be 
improved by using the approximate solution. 
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ECCENTRICITY 




5* (Continued) 

The analytical coast solution can be employed for three dimensional 
problems as well as planar cases. It is only necessary to rotate into 
the plane of motion at the start of coast and rotate back out of the 
plane to the original reference at the end of coast. The axes used must 
coincide with the perigee oriented axes used in developing the coast 
solution. 
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